#Background At Blue Oak Ranch Reserve, we established a 10m x 26m fenced field site to test whether evolution over the course of invasion away from roads has resulted in enhanced performance in undisturbed vegetation relative to roadside populations. The experiment was replicated in a randomized block design (20 plots in total). Each plot was 1.5m2 with 16 D. graveolens growing in a 4 x 4 grid centered on the plot. There was 33cm between each plant and 25cm between the edge plants and the border of the plot.
The experiment included multiple treatments; however, only the two most relevant to the focus of this paper are included here. We tested whether plant genotypes collected from the two habitats responded differently to the disturbance of biomass removal. We tilled in December 2020 to completely remove below and aboveground biomass, and then weeded to remove aboveground biomass throughout the growing season. In contrast, we left the control plots untouched, allowing the previous year’s thatch to persist and background vegetation to grow throughout the experiment.
In January 2021, we germinated seeds in Petri dishes at the UCSC Coastal Science Campus greenhouse incubation chambers before transplanting them into field-collected soil (collected in late December 2020 from Blue Oak Ranch Reserve). Seedlings grew in the greenhouse for about eight weeks until all plants had their first two true leaves emerge and lengthen. Ideally, we would have placed seeds directly into the field, but to maximize biosafety, we used seedling transplants that could be tracked with 100% certainty.
We measured the longest leaf for each plant and then transplanted them into the ground in late February 2021 at Blue Oak Ranch Reserve. During the first month of growth, we replaced any D. graveolens that died. We conducted weekly phenology surveys to assess D. graveolens plant health, and at the first sign of buds, we measured plant height and harvested the aboveground biomass by cutting at the root crown and drying in a 60ºC oven for 3 days before weighing.
#Data Analysis Statistical analyses were performed in R version 4.2.1 (R Core Team 2022) using linear mixed-effects models with the lme4 (Bates et al. 2015), lmerTest (Kuznetsova et al. 2017), and DHARMa packages (Hartig 2022), generalized linear mixed models with the glmmTMB package (Brooks et al. 2017), and mixed effects cox models with the coxme (Therneau 2022a) and survival (Therneau 2022b) packages.
##Models Included
###Cox proportional hazard models Here we will use ‘coxme’ which allows you to conduct mixed effects Cox proportional hazards models. Information here: https://cran.r-project.org/web/packages/coxme/vignettes/coxme.pdf. We conducted a germination experiment using Dittrichia graveolens seeds on filter paper. “Surv” creates a survival object to combine the days column (NumDaysAlive) and the censor column (Censor) to be used as a response variable in a model formula.The model follows the same syntax as linear models (lm, lmer, etc). Fixed effects: “var1 * var2” will give you the interaction term and the individual variables: so “var1 + var2 + var1:var2”. To add random effects, type “+ (1|random effect)”.
Assumptions for cox models: https://www.theanalysisfactor.com/assumptions-cox-regression/
###Linear mixed-effects models Linear Mixed Effects models are used for regression analyses involving dependent data. For a tutorial: https://doi.org/10.1177/2515245920960351
fullmodel<-lmer(log(Biomass)~ #Response variable: biomass HabitatTreatment+ #Fixed effects and their interactions() (1|Site)+(1|Block), #Random effect with random intercept only data=mydata) #Dataframe
###Generalized linear mixed models under construction
#Libraries
#install.packages("coxme")
#install.packages("survival")
#install.packages("ggplot2")
#install.packages("ggfortify")
#install.packages("car")
#install.packages("multcomp")
#install.packages("lme4")
#install.packages("lmerTest")
#install.packages("DHARMa")
#install.packages("dplyr")
#install.packages("emmeans")
#install.packages('TMB', type = 'source')
#install.packages("glmmTMB")
#install.packages("MASS")
#install.packages("emmeans")
#install.packages("AICcmodavg")
library(coxme)
library(survival)
library(ggplot2)
library(ggfortify)
library(car)
library(multcomp)
library(lme4)
library(lmerTest)
library(DHARMa)
library(dplyr)
library(emmeans)
library(TMB)
library(glmmTMB)
library(MASS)
library(emmeans)
library(AICcmodavg)
#Load Data This dataframe has one row per plant (800 observations). Data are for survivorship curves (3 censor options), the number of days the plant stayed alive (NumDaysAlive) and aboveground biomass. Censors with a 1 denote reaching the event (CensorAll = died, CensorBiomass = survived to collect biomass, CensorReproduction = survived to reproduce) and a 0 denoting when a seed didn’t germinate by the last census date (Census = 11/15/21). CensorReproduction will be most useful in understanding the amount of biomass produced by an individual when buds appear.
str(mydata) #Check that each column has the right class (factor, integer, numeric, etc.)
'data.frame': 800 obs. of 29 variables:
$ Block : Factor w/ 10 levels "A","B","C","D",..: 1 2 3 4 5 6 7 8 9 10 ...
$ Plot : int 5 2 2 1 5 1 1 4 4 2 ...
$ Flag : Factor w/ 50 levels "A1","A2","A3",..: 5 7 12 16 25 26 31 39 44 47 ...
$ Pos : int 3 15 13 13 6 10 8 4 12 3 ...
$ Flag_Pos : Factor w/ 800 levels "A1_01","A1_02",..: 67 111 189 253 390 410 488 612 700 739 ...
$ Population : Factor w/ 16 levels "BAY-A","BAY-O",..: 1 1 1 1 1 1 1 1 1 1 ...
$ Site : Factor w/ 8 levels "BAY","CHE","GUA",..: 1 1 1 1 1 1 1 1 1 1 ...
$ Habitat : Factor w/ 2 levels "Roadside","Vegetated": 1 1 1 1 1 1 1 1 1 1 ...
$ Treatment : Factor w/ 5 levels "Biomass Removal",..: 1 1 1 1 1 1 1 1 1 1 ...
$ PlantDate : Factor w/ 8 levels "2/27/21","3/1/21",..: 5 1 1 1 1 1 8 1 1 1 ...
$ LeafMeas1 : num 30 20.3 27 14.9 16.7 ...
$ LeafMeas2 : int 46 51 NA 56 42 45 36 58 52 45 ...
$ Growth : num 16 30.7 -27 41 25.3 ...
$ MortDate : Factor w/ 29 levels "","10/1/21","11/14/21",..: 1 25 7 1 1 1 1 1 1 1 ...
$ MortHarvDate : Factor w/ 33 levels "10/1/21","10/28/21",..: 21 27 9 29 20 29 29 29 30 30 ...
$ CensorAll : int 1 1 1 1 1 1 1 1 1 1 ...
$ DaysMort : int 105 181 62 195 139 195 185 195 198 198 ...
$ Census : Factor w/ 1 level "11/15/21": 1 1 1 1 1 1 1 1 1 1 ...
$ DaysCensus : int 241 261 261 261 261 261 251 261 261 261 ...
$ NumDaysAlive : int 105 181 62 195 139 195 185 195 198 198 ...
$ HarvDate : Factor w/ 22 levels "","10/1/21","10/28/21",..: 8 15 1 17 7 17 17 17 18 18 ...
$ CensorBiomass : int 1 1 0 1 1 1 1 1 1 1 ...
$ BudDate : Factor w/ 20 levels "","10/1/21","10/28/21",..: 9 1 1 17 11 17 17 17 17 17 ...
$ CensorReproduction: int 1 0 0 1 1 1 1 1 1 1 ...
$ PropBud : num 0.763 0.763 0.763 0.763 0.763 0.763 0.763 0.763 0.763 0.763 ...
$ PropBudSite : num 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 ...
$ Height : num 12.4 33.8 NA 56.6 21.5 55.2 45.2 37.5 50.8 48.5 ...
$ Biomass : num 0.762 7.048 NA 7.803 2.572 ...
$ Biomass.date : Factor w/ 18 levels "","0.561","10/13/21",..: 17 16 1 8 15 11 11 11 10 10 ...
#Early Growth This code uses Growth data with Habitat (roadside and vegetated) and Treatment in a glmm model. Anova and Tukey tests are used on the successful Model4 with the creation of a box plot as a finished product.
Note: 10 blocks, 5 treatments, 16 populations (CHE-O), 8 sites (random: population pairs; CHE), 2 habitat (fixed effect: roadside and vegetated; R or O)
Also Note: Growth data is a measure of growth of longest leaf. Each plant was measured upon transplanting into the field in March and then again in June 2021. This plant has a juvenile stage of a basal rosette, and then it bolts and produced smaller cauline leaves. The goal was to capture early growth data for this plant before bolting occurs so that the measurements capture the basal rosette stage, however in some cases the plants bolted earlier than expected and the resulting measurement was smaller than previous. In these cases we determined that the data should be removed from the dataset as the negative number (or changing it to a zero) does not reflect the biological importance of the measurement.
Here we need to filter the data to remove Growth>0 and to convert plant measurement dates to date format
growth_mydata<-mydata%>%filter(Growth>0) #Here I am only looking at the Growth data that is greater than 0 (see Also Note above)
growth_mydata$PlantDate<-as.Date(growth_mydata$PlantDate,"%m/%d/%y")
growth_mydata$Num.Days.Growth<-as.Date("2021-05-22")-growth_mydata$PlantDate
growth_mydata$Num.Days.Growth<-as.numeric(growth_mydata$Num.Days.Growth)
growth_mydata$Growth.Rate<-growth_mydata$Growth/growth_mydata$Num.Days.Growth #First calculate number of days of growth to get the Rate (Growth/Num.Days.Growing). Then fit data to Beta distribution
##Histograms Original data
hist(growth_mydata$Growth.Rate,
col='steelblue',main='Original') #Original data is skewed, let's test for normality and consider log transforming the data
shapiro.test(growth_mydata$Growth.Rate)
Shapiro-Wilk normality test
data: growth_mydata$Growth.Rate
W = 0.90445, p-value < 2.2e-16
Log transform data (https://www.statology.org/transform-data-in-r/)
log_growth_mydata<-log10(growth_mydata$Growth.Rate)
hist(log_growth_mydata,
col='steelblue',main='Log Transformed') #Log transformed data, this looks better than the original distribution
shapiro.test(log_growth_mydata) #Data does not improve with log transformation
Shapiro-Wilk normality test
data: log_growth_mydata
W = 0.92178, p-value = 1.503e-15
Next I tried a square root transformation, which improved the distribution, but my first model failed the singular fit (see Model 1).
sqrt_growth_mydata<-sqrt(growth_mydata$Growth.Rate)
hist(sqrt_growth_mydata,
col='steelblue',main='Log Transformed') #Square root transformed data looks better than the original distribution
shapiro.test(sqrt_growth_mydata)
Shapiro-Wilk normality test
data: sqrt_growth_mydata
W = 0.98581, p-value = 7.695e-05
##Model 1 - lmer: Modeling Habitat and Treatment with Site and Block as random effects
growth_fullmodel1<-lmer(sqrt(Growth.Rate)~Habitat*Treatment+(1|Site)+(1|Block),data=growth_mydata)
boundary (singular) fit: see help('isSingular')
isSingular(growth_fullmodel1,tol=1e-4) #=True
[1] TRUE
summary(growth_fullmodel1) #Variance explained by Site = 0.000
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: sqrt(Growth.Rate) ~ Habitat * Treatment + (1 | Site) + (1 | Block)
Data: growth_mydata
REML criterion at convergence: -583.3
Scaled residuals:
Min 1Q Median 3Q Max
-2.61595 -0.69539 0.05978 0.73053 2.88562
Random effects:
Groups Name Variance Std.Dev.
Block (Intercept) 0.001598 0.03998
Site (Intercept) 0.000000 0.00000
Residual 0.016171 0.12716
Number of obs: 506, groups: Block, 10; Site, 8
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 0.271164 0.024328 78.070400 11.146 < 2e-16
HabitatVegetated 0.023733 0.028530 487.902365 0.832 0.40590
TreatmentBiomass Removal 0.299248 0.025549 488.988743 11.713 < 2e-16
TreatmentHemizonia 0.009046 0.029694 488.862262 0.305 0.76078
TreatmentRaking 0.037015 0.028063 487.757926 1.319 0.18778
TreatmentRaking & Clipping 0.073071 0.027082 488.008833 2.698 0.00722
HabitatVegetated:TreatmentBiomass Removal -0.005098 0.035439 487.461999 -0.144 0.88568
HabitatVegetated:TreatmentHemizonia 0.014342 0.041419 487.869509 0.346 0.72928
HabitatVegetated:TreatmentRaking -0.022588 0.038582 488.126281 -0.585 0.55852
HabitatVegetated:TreatmentRaking & Clipping -0.006018 0.037436 488.124609 -0.161 0.87236
(Intercept) ***
HabitatVegetated
TreatmentBiomass Removal ***
TreatmentHemizonia
TreatmentRaking
TreatmentRaking & Clipping **
HabitatVegetated:TreatmentBiomass Removal
HabitatVegetated:TreatmentHemizonia
HabitatVegetated:TreatmentRaking
HabitatVegetated:TreatmentRaking & Clipping
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) HbttVg TrtmBR TrtmnH TrtmnR TrtR&C HV:TBR HbV:TH HbV:TR
HabittVgttd -0.614
TrtmntBmssR -0.694 0.584
TretmntHmzn -0.596 0.501 0.567
TretmntRkng -0.626 0.532 0.597 0.512
TrtmntRkn&C -0.650 0.551 0.619 0.533 0.560
HbttVgt:TBR 0.493 -0.804 -0.714 -0.402 -0.428 -0.442
HbttVgtt:TH 0.418 -0.687 -0.398 -0.708 -0.365 -0.378 0.553
HbttVgtt:TR 0.454 -0.741 -0.432 -0.369 -0.726 -0.406 0.596 0.508
HbttVg:TR&C 0.465 -0.761 -0.442 -0.381 -0.403 -0.721 0.612 0.526 0.564
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
anova(growth_fullmodel1)
Type III Analysis of Variance Table with Satterthwaite's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
Habitat 0.0464 0.04642 1 488.79 2.8709 0.09083 .
Treatment 7.4582 1.86456 4 489.79 115.3056 < 2e-16 ***
Habitat:Treatment 0.0146 0.00366 4 487.95 0.2260 0.92382
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Site as a random effect does not explain any of the variance in the model, therefore let’s try Site as a fixed effect to see if it adds to the model.
##Model 2 - lmer: Modeling Site as a fixed effect and Block as a random effect
growth_fullmodel2<-lmer(log(Growth.Rate)~
Habitat*Treatment+
Site+
(1|Block),
data=growth_mydata)
summary(growth_fullmodel2) #As a fixed effect, one of the Sites (OAP) is significant.
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: log(Growth.Rate) ~ Habitat * Treatment + Site + (1 | Block)
Data: growth_mydata
REML criterion at convergence: 1374.5
Scaled residuals:
Min 1Q Median 3Q Max
-6.2000 -0.4782 0.1779 0.6492 1.8316
Random effects:
Groups Name Variance Std.Dev.
Block (Intercept) 0.06789 0.2606
Residual 0.82580 0.9087
Number of obs: 506, groups: Block, 10
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) -2.984938 0.204110 165.530302 -14.624 < 2e-16
HabitatVegetated 0.220653 0.204078 480.955621 1.081 0.28014
TreatmentBiomass Removal 1.729151 0.183091 482.054827 9.444 < 2e-16
TreatmentHemizonia 0.196272 0.212666 481.942503 0.923 0.35651
TreatmentRaking 0.346745 0.201477 480.773125 1.721 0.08589
TreatmentRaking & Clipping 0.547314 0.194472 480.894171 2.814 0.00509
SiteCHE 0.039617 0.163190 481.278240 0.243 0.80829
SiteGUA -0.006678 0.163292 481.229118 -0.041 0.96740
SiteLEX -0.063730 0.164176 481.225590 -0.388 0.69805
SiteOAP 0.329379 0.161428 481.228603 2.040 0.04186
SitePAR 0.105935 0.165141 480.587844 0.641 0.52152
SitePEN 0.115262 0.167219 481.510246 0.689 0.49098
SiteSSJ 0.112886 0.166554 480.788469 0.678 0.49824
HabitatVegetated:TreatmentBiomass Removal -0.152103 0.253524 480.438799 -0.600 0.54882
HabitatVegetated:TreatmentHemizonia 0.048304 0.296429 480.886893 0.163 0.87062
HabitatVegetated:TreatmentRaking -0.214508 0.276049 481.225843 -0.777 0.43750
HabitatVegetated:TreatmentRaking & Clipping -0.193752 0.268159 481.243075 -0.723 0.47032
(Intercept) ***
HabitatVegetated
TreatmentBiomass Removal ***
TreatmentHemizonia
TreatmentRaking .
TreatmentRaking & Clipping **
SiteCHE
SiteGUA
SiteLEX
SiteOAP *
SitePAR
SitePEN
SiteSSJ
HabitatVegetated:TreatmentBiomass Removal
HabitatVegetated:TreatmentHemizonia
HabitatVegetated:TreatmentRaking
HabitatVegetated:TreatmentRaking & Clipping
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation matrix not shown by default, as p = 17 > 12.
Use print(x, correlation=TRUE) or
vcov(x) if you need it
anova(growth_fullmodel2) #Site as a fixed effect accounts for 3% of the variance
Type III Analysis of Variance Table with Satterthwaite's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
Habitat 1.636 1.636 1 482.01 1.9816 0.1599
Treatment 211.245 52.811 4 483.04 63.9513 <2e-16 ***
Site 6.623 0.946 7 481.06 1.1457 0.3331
Habitat:Treatment 1.190 0.297 4 481.03 0.3601 0.8370
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Now we’ll look at the QQ plots and the residuals using DHARMa
qqnorm(resid(growth_fullmodel2)) #qqplot
qqline(resid(growth_fullmodel2)) #add the line
testDispersion(growth_fullmodel2) #red line should be in the middle of the distribution
DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated
data: simulationOutput
dispersion = 0.97037, p-value = 0.68
alternative hypothesis: two.sided
myDHARMagraph2<-simulateResiduals(growth_fullmodel2) #making a graph using DHARMa package, also testing for heteroscedasticity, https://cran.r-project.org/web/packages/DHARMa/vignettes/DHARMa.html#heteroscedasticity
plot(myDHARMagraph2) #plotting graph. At this point, you don't want any text or lines to be red. The QQ plots are not looking good for this model. Let's try fitting to a negative binomial distribution.
##Model 3 - glmer.nb: Negative Binomial
So this model isn’t working well either. Let’s try building a model and fitting it to a Beta distribution, first without a link function. Check it with DHARMa, and if it doesn’t look good, then fit it to a Beta distribution with a logit link function.
##Model 4 - glmm: Beta Distribution Now I’m using glmm because I’m fitting to other distributions (beta)
growth_fullmodel4<-glmmTMB(Growth.Rate~
Habitat*Treatment+
(1|Site)+
(1|Block),
family=beta_family(),
data=growth_mydata)
summary(growth_fullmodel4)
Family: beta ( logit )
Formula: Growth.Rate ~ Habitat * Treatment + (1 | Site) + (1 | Block)
Data: growth_mydata
AIC BIC logLik deviance df.resid
-989.3 -934.4 507.7 -1015.3 493
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
Site (Intercept) 0.0002645 0.01626
Block (Intercept) 0.0429336 0.20720
Number of obs: 506, groups: Site, 8; Block, 10
Dispersion parameter for beta family (): 11.4
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.27681 0.14929 -15.251 <2e-16 ***
HabitatVegetated 0.12961 0.17864 0.726 0.468
TreatmentBiomass Removal 1.60566 0.15070 10.655 <2e-16 ***
TreatmentHemizonia 0.07352 0.18737 0.392 0.695
TreatmentRaking 0.20030 0.17566 1.140 0.254
TreatmentRaking & Clipping 0.41358 0.16659 2.483 0.013 *
HabitatVegetated:TreatmentBiomass Removal -0.03513 0.20401 -0.172 0.863
HabitatVegetated:TreatmentHemizonia 0.10730 0.25788 0.416 0.677
HabitatVegetated:TreatmentRaking -0.09786 0.23867 -0.410 0.682
HabitatVegetated:TreatmentRaking & Clipping -0.06623 0.22709 -0.292 0.771
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
qqnorm(resid(growth_fullmodel4)) #qqplot
qqline(resid(growth_fullmodel4)) #add the line
testDispersion(growth_fullmodel4) #red line should be in the middle of the distribution
DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated
data: simulationOutput
dispersion = 0.97597, p-value = 0.84
alternative hypothesis: two.sided
myDHARMagraph4<-simulateResiduals(growth_fullmodel4) #making a graph using DHARMa package, also testing for heteroscedasticity, https://cran.r-project.org/web/packages/DHARMa/vignettes/DHARMa.html#heteroscedasticity
plot(myDHARMagraph4) #plotting graph. Looks good, but check the outliers to make sure they are real.
We should check the outliers in the the DHARMa plot to make sure they make sense.
growth_outlier_boxplot1<-ggplot(growth_mydata)+
geom_boxplot(aes(x=Habitat,
y=Growth.Rate),
size=0.5)+
theme_classic()+
ylim(0,1)+
scale_fill_manual(values=c("forestgreen","gray85"))+
labs(y="Early Growth Rate",
fill="Habitat",
x="Habitat")
growth_outlier_boxplot1
max(growth_mydata$Growth.Rate[growth_mydata$Habitat=="Away"])
Warning: no non-missing arguments to max; returning -Inf
[1] -Inf
max(growth_mydata$Growth.Rate[growth_mydata$Habitat=="Roadside"]) #Yes, these outliers make sense, so I don't need to worry about the red stars in the DHARMa plot.
[1] 0.75
###Best Model
growth_fullmodel4<-glmmTMB(Growth.Rate~
Habitat*Treatment+
(1|Site)+
(1|Block),
family=beta_family(),
data=growth_mydata)
summary(growth_fullmodel4)
Family: beta ( logit )
Formula: Growth.Rate ~ Habitat * Treatment + (1 | Site) + (1 | Block)
Data: growth_mydata
AIC BIC logLik deviance df.resid
-989.3 -934.4 507.7 -1015.3 493
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
Site (Intercept) 0.0002645 0.01626
Block (Intercept) 0.0429336 0.20720
Number of obs: 506, groups: Site, 8; Block, 10
Dispersion parameter for beta family (): 11.4
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.27681 0.14929 -15.251 <2e-16 ***
HabitatVegetated 0.12961 0.17864 0.726 0.468
TreatmentBiomass Removal 1.60566 0.15070 10.655 <2e-16 ***
TreatmentHemizonia 0.07352 0.18737 0.392 0.695
TreatmentRaking 0.20030 0.17566 1.140 0.254
TreatmentRaking & Clipping 0.41358 0.16659 2.483 0.013 *
HabitatVegetated:TreatmentBiomass Removal -0.03513 0.20401 -0.172 0.863
HabitatVegetated:TreatmentHemizonia 0.10730 0.25788 0.416 0.677
HabitatVegetated:TreatmentRaking -0.09786 0.23867 -0.410 0.682
HabitatVegetated:TreatmentRaking & Clipping -0.06623 0.22709 -0.292 0.771
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
###Post-Hoc Test Remove non-significant interaction terms before running the Tukey.
#?emmeans,emmeans(model,pairwise~treatment)
growth_fullmodel4.1<-glmmTMB(Growth.Rate~
Treatment+
(1|Site)+
(1|Block),
family=beta_family(),
data=growth_mydata)
summary(growth_fullmodel4.1)
Family: beta ( logit )
Formula: Growth.Rate ~ Treatment + (1 | Site) + (1 | Block)
Data: growth_mydata
AIC BIC logLik deviance df.resid
-996.0 -962.2 506.0 -1012.0 498
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
Site (Intercept) 0.0006051 0.0246
Block (Intercept) 0.0416553 0.2041
Number of obs: 506, groups: Site, 8; Block, 10
Dispersion parameter for beta family (): 11.3
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.2066 0.1144 -19.284 < 2e-16 ***
TreatmentBiomass Removal 1.5835 0.1053 15.043 < 2e-16 ***
TreatmentHemizonia 0.1210 0.1295 0.934 0.35019
TreatmentRaking 0.1507 0.1203 1.253 0.21035
TreatmentRaking & Clipping 0.3781 0.1150 3.287 0.00101 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
emmeans(growth_fullmodel4.1,
pairwise~Treatment)
$emmeans
Treatment emmean SE df lower.CL upper.CL
Grassland -2.207 0.1144 498 -2.431 -1.982
Biomass Removal -0.623 0.0818 498 -0.784 -0.462
Hemizonia -2.086 0.1158 498 -2.313 -1.858
Raking -2.056 0.1050 498 -2.262 -1.850
Raking & Clipping -1.828 0.0969 498 -2.019 -1.638
Results are given on the logit (not the response) scale.
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Grassland - Biomass Removal -1.5835 0.1053 498 -15.043 <.0001
Grassland - Hemizonia -0.1210 0.1295 498 -0.934 0.8835
Grassland - Raking -0.1507 0.1203 498 -1.253 0.7203
Grassland - Raking & Clipping -0.3781 0.1150 498 -3.287 0.0095
Biomass Removal - Hemizonia 1.4625 0.1071 498 13.658 <.0001
Biomass Removal - Raking 1.4327 0.0950 498 15.081 <.0001
Biomass Removal - Raking & Clipping 1.2054 0.0864 498 13.956 <.0001
Hemizonia - Raking -0.0298 0.1225 498 -0.243 0.9992
Hemizonia - Raking & Clipping -0.2571 0.1170 498 -2.198 0.1822
Raking - Raking & Clipping -0.2274 0.1062 498 -2.141 0.2043
Results are given on the log odds ratio (not the response) scale.
P value adjustment: tukey method for comparing a family of 5 estimates
###ggplot ####Interaction plot #####Biomass x Treatment: Growth
pd<-position_dodge(0)
growth_mydata1<-growth_mydata%>%
replace_na(list(Biomass=0))%>%
filter(CensorReproduction==1)%>%
filter(Treatment=="Biomass Removal"|Treatment=="Grassland")%>%
group_by(Treatment,Habitat)%>%
summarise(Mean=mean(Biomass),
SD=sd(Biomass),
N=length(Biomass))%>%
mutate(SE=SD/sqrt(N))
`summarise()` has grouped output by 'Treatment'. You can override using the `.groups` argument.
field_int_bio_grow_gg<-ggplot(growth_mydata1,
aes(x=factor(Treatment,levels=c("Biomass Removal","Grassland")),
y=Mean,
shape=Habitat,
group=Habitat,
color=Habitat))+
ylab("Biomass (g)")+
coord_cartesian(ylim = c(0,10))+
xlab("")+
#ggtitle("Biomass of Dittrichia graveolens that budded")+
theme_bw(base_size=16)+
theme(panel.border=element_blank(),
panel.grid.major=element_blank(),
panel.grid.minor=element_blank(),
axis.line=element_line(colour="black"),
axis.title.y=ggtext::element_markdown(),
legend.position = c(0.8, 0.8))+
geom_line(size=1,
position=pd)+
geom_point(size=5,
position=pd)+
scale_color_manual(values=c("gray85","forestgreen"))+
geom_errorbar(width=0.15,
position=pd,
aes(ymin=(Mean-SE),
ymax=(Mean+SE)))+
labs(y="*Dittrichia graveolens* Biomass (g)",
color="Source Habitat",
shape="Source Habitat")
field_int_bio_grow_gg
ggsave(plot=field_int_bio_grow_gg,file="/Users/Miranda/Documents/Education/UC Santa Cruz/Dittrichia/Dittrichia_Analysis/figures/field_int_bio_grow_gg.png",width=25,height=15,units="cm",dpi=800)
#####Growth Rate x Treatment: Growth
pd<-position_dodge(0)
growth_mydata2<-growth_mydata%>%
replace_na(list(Growth.Rate=0))%>%
filter(CensorReproduction==1)%>%
filter(Treatment=="Biomass Removal"|Treatment=="Grassland")%>%
group_by(Treatment,Habitat)%>%
summarise(Mean=mean(Growth.Rate),
SD=sd(Growth.Rate),
N=length(Growth.Rate))%>%
mutate(SE=SD/sqrt(N))
`summarise()` has grouped output by 'Treatment'. You can override using the `.groups` argument.
field_int_rate_grow_gg<-ggplot(growth_mydata2,
aes(x=factor(Treatment,levels=c("Biomass Removal","Grassland")),
y=Mean,
shape=Habitat,
group=Habitat,
color=Habitat))+
coord_cartesian(ylim = c(0,0.5))+
xlab("")+
theme_bw(base_size=16)+
theme(panel.border=element_blank(),
panel.grid.major=element_blank(),
panel.grid.minor=element_blank(),
axis.line=element_line(colour="black"),
axis.title.y=ggtext::element_markdown(),
legend.position = c(0.8, 0.8))+
# ggtitle("Growth of Dittrichia graveolens that budded")+
ylab("*Dittrichia graveolens* \n (Change in Leaf Length/Days)")+
geom_line(size=1,
position=pd)+
geom_point(size=5,
position=pd)+
scale_color_manual(values=c("gray85","forestgreen"))+
geom_errorbar(width=0.15,
position=pd,
aes(ymin=(Mean-SE),
ymax=(Mean+SE)))+
labs(y="Change in *Dittrichia graveolens* Leaf Length/Day",
color="Source Habitat",
shape="Source Habitat")
field_int_rate_grow_gg
ggsave(plot=field_int_rate_grow_gg,file="/Users/Miranda/Documents/Education/UC Santa Cruz/Dittrichia/Dittrichia_Analysis/figures/field_int_rate_grow_gg.png",width=25,height=15,units="cm",dpi=800)
####Boxplot
#Reproductive Biomass This code uses Biomass data with Habitat (roadside and vegetated) and Treatment in a lmer model. ANOVA and Tukey are used on the successful Model 3 with the creation of a box plot as a finished product.
Note: 10 blocks, 5 treatments, 16 populations (CHE-O), 8 sites (random: population pairs; CHE), 2 habitat (fixed effect: roadside and vegetated; R or O)
Also Note: Biomass data is a measure of aboveground biomass of all individuals harvested from the site. In some cases this was before they budded, but we harvested the wilted plant in case that information was needed in the future. Here we only want to look at biomass (well, log(Biomass)) for reproductive individuals.
Here we need to subset the data to only look at Biomass when CensorReproduction = 1, so that all Biomass data is for reproductive individuals only.
reproduction_mydata<-subset(mydata,CensorReproduction%in%c('1')) #Here I am only looking at the Biomass data where the CensorReproduction = 1
##Histograms
hist(reproduction_mydata$Biomass,col='steelblue',main='Original') #Original data is skewed, let's test for normality and consider log transforming the data
shapiro.test(reproduction_mydata$Biomass)
Shapiro-Wilk normality test
data: reproduction_mydata$Biomass
W = 0.62548, p-value < 2.2e-16
Log transform data (https://www.statology.org/transform-data-in-r/)
log_reproduction_mydata<-log10(reproduction_mydata$Biomass)
hist(log_reproduction_mydata,col='steelblue',main='Log Transformed') #Log transformed data, this looks better than the original distribution
shapiro.test(log_reproduction_mydata)
Shapiro-Wilk normality test
data: log_reproduction_mydata
W = 0.9648, p-value = 9.101e-06
#Log transformed data has a better distribution than the original data so we will use the log transformed data with our models
##Model 1 - lmer: Modeling Habitat and Treatment with Site and Block as random effects
fullmodel1<-lmer(log(Biomass)~
Habitat*Treatment+
(1|Site)+
(1|Block),
data=reproduction_mydata)
isSingular(fullmodel1,tol=1e-4)
[1] FALSE
summary(fullmodel1) #Variance explained by Site = 0.000
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: log(Biomass) ~ Habitat * Treatment + (1 | Site) + (1 | Block)
Data: reproduction_mydata
REML criterion at convergence: 675.3
Scaled residuals:
Min 1Q Median 3Q Max
-3.2251 -0.6836 0.0555 0.6095 2.9640
Random effects:
Groups Name Variance Std.Dev.
Block (Intercept) 0.1836330 0.42852
Site (Intercept) 0.0003818 0.01954
Residual 0.8322599 0.91228
Number of obs: 247, groups: Block, 10; Site, 8
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) -1.47119 0.25189 65.18355 -5.841 1.8e-07
HabitatVegetated 0.52617 0.31275 229.17635 1.682 0.0939
TreatmentBiomass Removal 3.16732 0.24296 227.53081 13.036 < 2e-16
TreatmentHemizonia -0.05905 0.31382 229.47907 -0.188 0.8509
TreatmentRaking 0.25358 0.41060 211.72838 0.618 0.5375
TreatmentRaking & Clipping 0.54416 0.31428 229.25370 1.731 0.0847
HabitatVegetated:TreatmentBiomass Removal -0.59897 0.35315 229.14685 -1.696 0.0912
HabitatVegetated:TreatmentHemizonia -0.30666 0.44327 228.66064 -0.692 0.4898
HabitatVegetated:TreatmentRaking -0.20197 0.55618 220.56602 -0.363 0.7168
HabitatVegetated:TreatmentRaking & Clipping -0.79028 0.43205 229.37117 -1.829 0.0687
(Intercept) ***
HabitatVegetated .
TreatmentBiomass Removal ***
TreatmentHemizonia
TreatmentRaking
TreatmentRaking & Clipping .
HabitatVegetated:TreatmentBiomass Removal .
HabitatVegetated:TreatmentHemizonia
HabitatVegetated:TreatmentRaking
HabitatVegetated:TreatmentRaking & Clipping .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) HbttVg TrtmBR TrtmnH TrtmnR TrtR&C HV:TBR HbV:TH HbV:TR
HabittVgttd -0.564
TrtmntBmssR -0.738 0.582
TretmntHmzn -0.558 0.439 0.582
TretmntRkng -0.445 0.351 0.463 0.342
TrtmntRkn&C -0.570 0.455 0.591 0.450 0.361
HbttVgt:TBR 0.500 -0.883 -0.675 -0.390 -0.309 -0.403
HbttVgtt:TH 0.392 -0.702 -0.404 -0.690 -0.246 -0.316 0.619
HbttVgtt:TR 0.319 -0.556 -0.331 -0.252 -0.723 -0.258 0.490 0.391
HbttVg:TR&C 0.413 -0.727 -0.425 -0.316 -0.262 -0.725 0.641 0.507 0.405
anova(fullmodel1)
Type III Analysis of Variance Table with Satterthwaite's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
Habitat 0.84 0.843 1 229.56 1.0133 0.3152
Treatment 483.42 120.855 4 228.83 145.2134 <2e-16 ***
Habitat:Treatment 3.83 0.957 4 226.65 1.1497 0.3340
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Site as a random effect does not explain any of the variance in the model, therefore let's try Site as a fixed effect to demonstrate that it doesn't add to the model.
predict(fullmodel1)
1 4 5 6 7 8 9 10
0.9783831 2.1411337 1.3028518 1.8487338 1.6247945 1.9815684 1.5191378 1.4884911
11 13 14 15 16 18 19 21
0.9832969 1.7787921 2.1460476 1.3077657 1.8536476 1.9864822 1.5240516 0.9837652
22 23 24 25 26 27 28 29
2.2719098 1.7792604 2.1465159 1.3082340 1.8541159 1.6301767 1.9869505 1.5245199
30 31 32 33 35 37 38 40
1.4938733 0.9834082 2.2715527 1.7789034 1.3078770 1.6298196 1.9865935 1.4935163
41 42 44 46 47 48 49 51
0.9827552 2.2708997 2.1455059 1.8531059 1.6291666 1.9859405 1.5235099 0.9783656
52 53 55 56 57 58 59 60
2.2665102 1.7738608 1.3028344 1.8487163 1.6247771 1.9815509 1.5191203 1.4884737
62 63 64 66 67 68 69 71
2.2710672 1.7784179 2.1456733 1.8532734 1.6293341 1.9861080 1.5236774 0.9827714
72 73 76 77 78 81 83 85
2.2709159 1.7782665 1.8531221 1.6291828 1.9859567 0.9055791 1.7010743 1.2300479
86 87 88 89 90 93 94 95
1.7759298 1.5519905 1.9087644 1.4463338 1.4156872 1.7059881 2.0732436 1.2349617
96 97 98 100 102 103 105 106
1.7808436 1.5569043 1.9136782 1.4206010 2.1991058 1.7064564 1.2354300 1.7813120
107 109 111 112 113 114 116 118
1.5573727 1.4517160 0.9106042 2.1987488 1.7060994 2.0733549 1.7809549 1.9137895
119 120 122 123 124 127 128 129
1.4513589 1.4207123 2.1980958 1.7054464 2.0727019 1.5563627 1.9131366 1.4507059
130 131 133 134 135 136 137 138
1.4200593 0.9055616 1.7010568 2.0683123 1.2300304 1.7759124 1.5519731 1.9087470
139 140 141 142 143 144 147 148
1.4463164 1.4156697 0.9101187 2.1982633 1.7056139 2.0728694 1.5565301 1.9133040
149 151 153 154 155 156 157 158
1.4508734 0.9099674 1.7054626 2.0727181 1.2344362 1.7803181 1.5563788 1.9131527
159 160 161 163 164 170 172 174
1.4507221 1.4200755 -2.1889361 -1.3934409 -1.0261854 -1.6788280 -0.8958777 -1.0212716
184 185 189 195 196 201 210 211
-1.0208033 -1.8590852 -1.6427992 -1.8594422 -1.3135603 -2.1845640 -1.6744559 -2.1889536
215 224 226 228 236 254 259 260
-1.8644848 -1.0216458 -1.3140458 -1.1812112 -1.3141971 -0.4951009 -1.1170969 -1.1477435
270 272 280 281 284 285 286 287
-1.1472751 -0.3695957 -1.1476322 -1.6583933 -0.4956426 -1.3339245 -0.7880426 -1.0119818
295 296 300 303 310 321 324 325
-1.3383141 -0.7924321 -1.1526748 -0.8627306 -1.1481177 -2.2479869 -1.0852362 -1.9235181
326 327 334 351 353 354 361 367
-1.3776362 -1.6015755 -1.0803224 -2.2429618 -1.4474666 -1.0802111 -2.2436148 -1.5972033
371 374 387 391 398 401 402 407
-2.2480044 -1.0852537 -1.5970359 -2.2435986 -1.2404133 -2.0284798 -0.7403353 -1.3820684
413 414 421 424 437 441 446 461
-1.2280708 -0.8608154 -2.0230977 -0.8603470 -1.3770433 -2.0241077 -1.1537570 -2.0239402
463 471 472 475 476 477 479 493
-1.2284451 -2.0240915 -0.7359470 -1.6996228 -1.1537408 -1.3776801 -1.4833368 -1.1349432
495 498 499 500 527 548 569 570
-1.6059696 -0.9272531 -1.3896837 -1.4203303 -1.2845687 -0.9276273 -1.0703998 -1.1010464
577 580 588 595 600 608 612 645
-0.9598293 -1.0961326 -0.6025870 -1.2816606 -1.0960213 -0.6035971 -0.3230274 -1.3203070
653 659 660 661 665 667 675 677
-0.8443668 -1.0991073 -1.1297539 -1.6393937 -1.3149249 -0.9929822 -1.3152819 -0.9933393
681 682 691 693 695 708 710 721
-1.6404037 -0.3522591 -1.6447933 -0.8492981 -1.3203245 -0.6370509 -1.1301281 -1.9088882
727 731 734 735 738 740 751 752
-1.2624768 -1.9039744 -0.7412237 -1.5795056 -0.9007890 -1.3938663 -1.9038631 -0.6157185
757 763 766 767 771 776 780 785
-1.2574516 -1.1090209 -1.0341654 -1.2581046 -1.9089056 -1.0385549 -1.3987976 -1.5798798
789 791 794 795 796 798 800
-1.3635939 -1.9044999 -0.7417492 -1.5800311 -1.0341492 -0.9013146 -1.3943918
hist(predict(fullmodel1,type="response"))
##Model 2 - lmer: Modeling Site as a fixed effect
fullmodel2<-lmer(log(Biomass)~
Habitat*Treatment+
Site+
(1|Block),
data=reproduction_mydata)
summary(fullmodel2)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: log(Biomass) ~ Habitat * Treatment + Site + (1 | Block)
Data: reproduction_mydata
REML criterion at convergence: 678.3
Scaled residuals:
Min 1Q Median 3Q Max
-3.2839 -0.6210 0.0399 0.5639 2.7969
Random effects:
Groups Name Variance Std.Dev.
Block (Intercept) 0.1650 0.4062
Residual 0.8361 0.9144
Number of obs: 247, groups: Block, 10
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) -1.720266 0.285577 105.107873 -6.024 2.54e-08
HabitatVegetated 0.488869 0.318875 223.009734 1.533 0.1267
TreatmentBiomass Removal 3.161011 0.244374 223.621849 12.935 < 2e-16
TreatmentHemizonia -0.004785 0.317523 223.398880 -0.015 0.9880
TreatmentRaking 0.156859 0.426831 224.250434 0.367 0.7136
TreatmentRaking & Clipping 0.509433 0.317070 223.835290 1.607 0.1095
SiteCHE 0.355322 0.232174 223.497740 1.530 0.1273
SiteGUA 0.410402 0.245081 221.592351 1.675 0.0954
SiteLEX 0.356617 0.233876 222.425736 1.525 0.1287
SiteOAP 0.316128 0.234233 222.925291 1.350 0.1785
SitePAR -0.001442 0.231840 221.753218 -0.006 0.9950
SitePEN 0.347827 0.244024 223.399330 1.425 0.1554
SiteSSJ 0.315678 0.237514 222.147337 1.329 0.1852
HabitatVegetated:TreatmentBiomass Removal -0.556028 0.358662 222.367673 -1.550 0.1225
HabitatVegetated:TreatmentHemizonia -0.363654 0.452306 222.295756 -0.804 0.4223
HabitatVegetated:TreatmentRaking -0.059139 0.573956 221.764455 -0.103 0.9180
HabitatVegetated:TreatmentRaking & Clipping -0.744475 0.440856 223.170989 -1.689 0.0927
(Intercept) ***
HabitatVegetated
TreatmentBiomass Removal ***
TreatmentHemizonia
TreatmentRaking
TreatmentRaking & Clipping
SiteCHE
SiteGUA .
SiteLEX
SiteOAP
SitePAR
SitePEN
SiteSSJ
HabitatVegetated:TreatmentBiomass Removal
HabitatVegetated:TreatmentHemizonia
HabitatVegetated:TreatmentRaking
HabitatVegetated:TreatmentRaking & Clipping .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation matrix not shown by default, as p = 17 > 12.
Use print(x, correlation=TRUE) or
vcov(x) if you need it
anova(fullmodel2)
Type III Analysis of Variance Table with Satterthwaite's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
Habitat 0.80 0.799 1 222.81 0.9554 0.3294
Treatment 478.33 119.583 4 224.94 143.0199 <2e-16 ***
Site 5.74 0.820 7 222.78 0.9811 0.4458
Habitat:Treatment 3.43 0.857 4 222.42 1.0244 0.3955
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Site as a fixed effect is not significant, therefore it should not be used as a fixed effect in addition to it not being used as a random effect.
##Model 3 - lmer: Site is removed from this model because it explains very little of the variance
fullmodel3<-lmer(log(Biomass)~
Habitat*Treatment+
(1|Block),
data=reproduction_mydata)
summary(fullmodel3)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: log(Biomass) ~ Habitat * Treatment + (1 | Block)
Data: reproduction_mydata
REML criterion at convergence: 675.3
Scaled residuals:
Min 1Q Median 3Q Max
-3.2235 -0.6825 0.0572 0.6111 2.9658
Random effects:
Groups Name Variance Std.Dev.
Block (Intercept) 0.1839 0.4288
Residual 0.8326 0.9124
Number of obs: 247, groups: Block, 10
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) -1.47141 0.25186 65.20283 -5.842 1.79e-07
HabitatVegetated 0.52673 0.31273 229.69739 1.684 0.0935
TreatmentBiomass Removal 3.16743 0.24299 230.30671 13.035 < 2e-16
TreatmentHemizonia -0.05973 0.31385 230.17407 -0.190 0.8492
TreatmentRaking 0.25499 0.41045 230.89163 0.621 0.5351
TreatmentRaking & Clipping 0.54463 0.31430 230.49078 1.733 0.0845
HabitatVegetated:TreatmentBiomass Removal -0.59959 0.35315 229.14496 -1.698 0.0909
HabitatVegetated:TreatmentHemizonia -0.30595 0.44325 229.12563 -0.690 0.4907
HabitatVegetated:TreatmentRaking -0.20398 0.55605 228.65197 -0.367 0.7141
HabitatVegetated:TreatmentRaking & Clipping -0.79088 0.43203 229.90026 -1.831 0.0685
(Intercept) ***
HabitatVegetated .
TreatmentBiomass Removal ***
TreatmentHemizonia
TreatmentRaking
TreatmentRaking & Clipping .
HabitatVegetated:TreatmentBiomass Removal .
HabitatVegetated:TreatmentHemizonia
HabitatVegetated:TreatmentRaking
HabitatVegetated:TreatmentRaking & Clipping .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) HbttVg TrtmBR TrtmnH TrtmnR TrtR&C HV:TBR HbV:TH HbV:TR
HabittVgttd -0.564
TrtmntBmssR -0.738 0.582
TretmntHmzn -0.558 0.439 0.582
TretmntRkng -0.445 0.351 0.463 0.342
TrtmntRkn&C -0.570 0.455 0.591 0.450 0.361
HbttVgt:TBR 0.500 -0.883 -0.675 -0.390 -0.309 -0.403
HbttVgtt:TH 0.392 -0.702 -0.405 -0.690 -0.246 -0.316 0.619
HbttVgtt:TR 0.319 -0.555 -0.331 -0.252 -0.723 -0.258 0.490 0.391
HbttVg:TR&C 0.413 -0.727 -0.425 -0.316 -0.262 -0.725 0.641 0.507 0.405
anova(fullmodel3)
Type III Analysis of Variance Table with Satterthwaite's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
Habitat 0.84 0.844 1 229.78 1.0140 0.3150
Treatment 483.49 120.873 4 231.54 145.1817 <2e-16 ***
Habitat:Treatment 3.83 0.959 4 229.30 1.1514 0.3332
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
###QQ Plots Now we’ll look at the QQ plots and the residuals using DHARMa
qqnorm(resid(fullmodel3)) #qqplot
qqline(resid(fullmodel3)) #add the line
testDispersion(fullmodel3) #red line should be in the middle of the distribution
DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated
data: simulationOutput
dispersion = 0.9764, p-value = 0.904
alternative hypothesis: two.sided
myDHARMagraph3<-simulateResiduals(fullmodel3) #making a graph using DHARMa package, also testing for heteroscedasticity, https://cran.r-project.org/web/packages/DHARMa/vignettes/DHARMa.html#heteroscedasticity
plot(myDHARMagraph3) #plotting graph. At this point, you don't want any text or lines to be red.
When we compare the model summaries and Anova of Models 2 & 3, we see that the removal of Site doesn’t impact the model. Therefore the more simple Model 3 is a better choice. But let’s keep checking…
###Compare AIC Scores Now I need to compare the AIC scores for all the models to tell me which is the better model (but it will not say which fits my data better, that is why I did all the DHARMa stuff)
Using aictab to make the comparison of models and table
models<-list(fullmodel1,fullmodel2,fullmodel3)
mod.names<-c('Site.Random','Site.Fixed','No.Site')
aictab(cand.set=models,modnames=mod.names)
Warning:
Model selection for fixed effects is only appropriate with ML estimation:
REML (default) should only be used to select random effects for a constant set of fixed effects
Model selection based on AICc:
K AICc Delta_AICc AICcWt Cum.Wt Res.LL
No.Site 12 700.60 0.00 0.75 0.75 -337.63
Site.Random 13 702.82 2.23 0.25 1.00 -337.63
Site.Fixed 19 719.64 19.04 0.00 1.00 -339.15
#The lowest AICc score is listed first and indicates the best fitting model, here, Model 3 (No.Site) where Site is not included. The cut-off for comparing Models is 2 units. The difference between Model 1 (Site.Random) and Model 3 (No.Site) is 2.12. So Model 3 is marginally better than Model 2.
##Other Models Attempted Other ideas that were explored to keep Site but resulted in singular error and Site explaining 0.000 of the variance: fullmodel4<-lmer(log(Biomass)~ Treatment + (1|Site) + (1|Block), data = mydata) summary(fullmodel4) fullmodel5<-lmer(log(Biomass)~ Treatment + Habitat + (1|Site) + (1|Block), data = mydata) summary(fullmodel5)
###Best Model
fullmodel3<-lmer(log(Biomass)~
Habitat*Treatment+
(1|Block),
data=reproduction_mydata)
summary(fullmodel3)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: log(Biomass) ~ Habitat * Treatment + (1 | Block)
Data: reproduction_mydata
REML criterion at convergence: 675.3
Scaled residuals:
Min 1Q Median 3Q Max
-3.2235 -0.6825 0.0572 0.6111 2.9658
Random effects:
Groups Name Variance Std.Dev.
Block (Intercept) 0.1839 0.4288
Residual 0.8326 0.9124
Number of obs: 247, groups: Block, 10
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) -1.47141 0.25186 65.20283 -5.842 1.79e-07
HabitatVegetated 0.52673 0.31273 229.69739 1.684 0.0935
TreatmentBiomass Removal 3.16743 0.24299 230.30671 13.035 < 2e-16
TreatmentHemizonia -0.05973 0.31385 230.17407 -0.190 0.8492
TreatmentRaking 0.25499 0.41045 230.89163 0.621 0.5351
TreatmentRaking & Clipping 0.54463 0.31430 230.49078 1.733 0.0845
HabitatVegetated:TreatmentBiomass Removal -0.59959 0.35315 229.14496 -1.698 0.0909
HabitatVegetated:TreatmentHemizonia -0.30595 0.44325 229.12563 -0.690 0.4907
HabitatVegetated:TreatmentRaking -0.20398 0.55605 228.65197 -0.367 0.7141
HabitatVegetated:TreatmentRaking & Clipping -0.79088 0.43203 229.90026 -1.831 0.0685
(Intercept) ***
HabitatVegetated .
TreatmentBiomass Removal ***
TreatmentHemizonia
TreatmentRaking
TreatmentRaking & Clipping .
HabitatVegetated:TreatmentBiomass Removal .
HabitatVegetated:TreatmentHemizonia
HabitatVegetated:TreatmentRaking
HabitatVegetated:TreatmentRaking & Clipping .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) HbttVg TrtmBR TrtmnH TrtmnR TrtR&C HV:TBR HbV:TH HbV:TR
HabittVgttd -0.564
TrtmntBmssR -0.738 0.582
TretmntHmzn -0.558 0.439 0.582
TretmntRkng -0.445 0.351 0.463 0.342
TrtmntRkn&C -0.570 0.455 0.591 0.450 0.361
HbttVgt:TBR 0.500 -0.883 -0.675 -0.390 -0.309 -0.403
HbttVgtt:TH 0.392 -0.702 -0.405 -0.690 -0.246 -0.316 0.619
HbttVgtt:TR 0.319 -0.555 -0.331 -0.252 -0.723 -0.258 0.490 0.391
HbttVg:TR&C 0.413 -0.727 -0.425 -0.316 -0.262 -0.725 0.641 0.507 0.405
###Post-Hoc Test You should remove non-significant interaction terms before running a post-hoc test. It is difficult to judge the main effects (Habitat and Site) when you also have the interaction term present, so when it is not significant, fit a new model without it.
fullmodel3.1<-lmer(log(Biomass)~
Treatment+
(1|Block),
data=reproduction_mydata)
summary(fullmodel3.1)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: log(Biomass) ~ Treatment + (1 | Block)
Data: reproduction_mydata
REML criterion at convergence: 677
Scaled residuals:
Min 1Q Median 3Q Max
-3.1914 -0.6482 0.0311 0.6048 3.1267
Random effects:
Groups Name Variance Std.Dev.
Block (Intercept) 0.1819 0.4265
Residual 0.8320 0.9121
Number of obs: 247, groups: Block, 10
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) -1.2310 0.2074 33.6195 -5.936 1.09e-06 ***
TreatmentBiomass Removal 2.8917 0.1790 236.6910 16.159 < 2e-16 ***
TreatmentHemizonia -0.1744 0.2261 237.1867 -0.772 0.441
TreatmentRaking 0.1873 0.2822 237.0309 0.664 0.508
TreatmentRaking & Clipping 0.1447 0.2140 235.2223 0.676 0.500
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) TrtmBR TrtmnH TrtmnR
TrtmntBmssR -0.673
TretmntHmzn -0.527 0.616
TretmntRkng -0.428 0.500 0.368
TrtmntRkn&C -0.553 0.644 0.516 0.412
anova(fullmodel3.1)
Type III Analysis of Variance Table with Satterthwaite's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
Treatment 490.29 122.57 4 236.51 147.33 < 2.2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
qqnorm(resid(fullmodel3.1)) #qqplot
qqline(resid(fullmodel3.1)) #add the line
testDispersion(fullmodel3.1) #red line should be in the middle of the distribution
DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated
data: simulationOutput
dispersion = 0.99335, p-value = 0.896
alternative hypothesis: two.sided
myDHARMagraph3.1<-simulateResiduals(fullmodel3.1) #testing for heteroscedasticity
plot(myDHARMagraph3.1) #plotting graph
fullmodel3.2<-lmer(log(Biomass)~
Treatment+
(1|Block)+
(1|Population),
data=mydata)
summary(fullmodel3.2)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: log(Biomass) ~ Treatment + (1 | Block) + (1 | Population)
Data: mydata
REML criterion at convergence: 1381.1
Scaled residuals:
Min 1Q Median 3Q Max
-3.4116 -0.6367 -0.0001 0.6389 3.6745
Random effects:
Groups Name Variance Std.Dev.
Population (Intercept) 0.000801 0.0283
Block (Intercept) 0.060607 0.2462
Residual 1.290867 1.1362
Number of obs: 441, groups: Population, 16; Block, 10
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) -2.0644 0.1565 52.0231 -13.195 <2e-16 ***
TreatmentBiomass Removal 3.5659 0.1650 426.7907 21.609 <2e-16 ***
TreatmentHemizonia -0.0984 0.1952 432.3439 -0.504 0.6145
TreatmentRaking -0.2757 0.1920 428.3744 -1.436 0.1517
TreatmentRaking & Clipping 0.3912 0.1822 428.2538 2.147 0.0324 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) TrtmBR TrtmnH TrtmnR
TrtmntBmssR -0.712
TretmntHmzn -0.602 0.571
TretmntRkng -0.608 0.577 0.484
TrtmntRkn&C -0.641 0.607 0.519 0.519
anova(fullmodel3.2)
Type III Analysis of Variance Table with Satterthwaite's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
Treatment 1233 308.25 4 426.32 238.79 < 2.2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
qqnorm(resid(fullmodel3.2)) #qqplot
qqline(resid(fullmodel3.2)) #add the line
testDispersion(fullmodel3.2) #red line should be in the middle of the distribution
DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated
data: simulationOutput
dispersion = 0.98522, p-value = 0.904
alternative hypothesis: two.sided
myDHARMagraph3.2<-simulateResiduals(fullmodel3.2) #testing for heteroscedasticity
plot(myDHARMagraph3.2) #plotting graph
Let’s use emmeans for our Best Model (fullmodel3) and the two off-shoots.
#?emmeans, emmeans(model, pairwise ~ treatment)
emmeans(fullmodel3,pairwise~Treatment)
NOTE: Results may be misleading due to involvement in interactions
$emmeans
Treatment emmean SE df lower.CL upper.CL
Grassland -1.21 0.209 33.3 -1.63 -0.784
Biomass Removal 1.66 0.159 11.9 1.31 2.006
Hemizonia -1.42 0.213 35.1 -1.85 -0.989
Raking -1.06 0.272 80.4 -1.60 -0.514
Raking & Clipping -1.06 0.202 29.9 -1.47 -0.646
Results are averaged over the levels of: Habitat
Degrees-of-freedom method: kenward-roger
Results are given on the log (not the response) scale.
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Grassland - Biomass Removal -2.8676 0.180 232 -15.931 <.0001
Grassland - Hemizonia 0.2127 0.228 232 0.934 0.8833
Grassland - Raking -0.1530 0.285 232 -0.537 0.9834
Grassland - Raking & Clipping -0.1492 0.217 230 -0.687 0.9592
Biomass Removal - Hemizonia 3.0803 0.183 233 16.789 <.0001
Biomass Removal - Raking 2.7146 0.249 231 10.886 <.0001
Biomass Removal - Raking & Clipping 2.7184 0.172 231 15.811 <.0001
Hemizonia - Raking -0.3657 0.292 234 -1.251 0.7211
Hemizonia - Raking & Clipping -0.3619 0.220 231 -1.645 0.4700
Raking - Raking & Clipping 0.0038 0.279 231 0.014 1.0000
Results are averaged over the levels of: Habitat
Degrees-of-freedom method: kenward-roger
Results are given on the log (not the response) scale.
P value adjustment: tukey method for comparing a family of 5 estimates
emmeans(fullmodel3.1,pairwise~Treatment)
$emmeans
Treatment emmean SE df lower.CL upper.CL
Grassland -1.23 0.208 33.4 -1.65 -0.809
Biomass Removal 1.66 0.158 12.0 1.32 2.006
Hemizonia -1.41 0.212 35.4 -1.84 -0.976
Raking -1.04 0.270 80.2 -1.58 -0.507
Raking & Clipping -1.09 0.199 29.0 -1.49 -0.678
Degrees-of-freedom method: kenward-roger
Results are given on the log (not the response) scale.
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Grassland - Biomass Removal -2.8917 0.179 237 -16.126 <.0001
Grassland - Hemizonia 0.1744 0.227 237 0.770 0.9391
Grassland - Raking -0.1873 0.283 237 -0.662 0.9642
Grassland - Raking & Clipping -0.1447 0.214 235 -0.675 0.9616
Biomass Removal - Hemizonia 3.0661 0.183 238 16.756 <.0001
Biomass Removal - Raking 2.7044 0.248 236 10.916 <.0001
Biomass Removal - Raking & Clipping 2.7469 0.169 236 16.249 <.0001
Hemizonia - Raking -0.3617 0.291 239 -1.245 0.7249
Hemizonia - Raking & Clipping -0.3192 0.217 236 -1.470 0.5829
Raking - Raking & Clipping 0.0425 0.276 236 0.154 0.9999
Degrees-of-freedom method: kenward-roger
Results are given on the log (not the response) scale.
P value adjustment: tukey method for comparing a family of 5 estimates
emmeans(fullmodel3.2,pairwise~Treatment)
$emmeans
Treatment emmean SE df lower.CL upper.CL
Grassland -2.06 0.157 57.5 -2.38 -1.75
Biomass Removal 1.50 0.122 23.1 1.25 1.75
Hemizonia -2.16 0.161 62.1 -2.49 -1.84
Raking -2.34 0.158 58.9 -2.66 -2.02
Raking & Clipping -1.67 0.146 44.0 -1.97 -1.38
Degrees-of-freedom method: kenward-roger
Results are given on the log (not the response) scale.
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Grassland - Biomass Removal -3.5659 0.166 428 -21.525 <.0001
Grassland - Hemizonia 0.0984 0.196 433 0.501 0.9873
Grassland - Raking 0.2757 0.193 429 1.429 0.6092
Grassland - Raking & Clipping -0.3912 0.183 429 -2.136 0.2067
Biomass Removal - Hemizonia 3.6643 0.170 428 21.602 <.0001
Biomass Removal - Raking 3.8416 0.166 426 23.076 <.0001
Biomass Removal - Raking & Clipping 3.1746 0.155 423 20.486 <.0001
Hemizonia - Raking 0.1773 0.198 433 0.896 0.8985
Hemizonia - Raking & Clipping -0.4896 0.186 427 -2.631 0.0666
Raking - Raking & Clipping -0.6669 0.184 423 -3.620 0.0030
Degrees-of-freedom method: kenward-roger
Results are given on the log (not the response) scale.
P value adjustment: tukey method for comparing a family of 5 estimates
#These models result in similar Tukey outcomes
##ggplot - lmer
ggplot(data=reproduction_mydata,
aes(x=Treatment,
y=log(Biomass)))+
geom_boxplot() #plot data from log(data)
ggplot(data=reproduction_mydata,
aes(x=Habitat,
y=log(Biomass)))+
geom_boxplot() #plot data from log(data)
#Survival Analysis This code uses NumDaysAlive data with Habitat (roadside and vegetated) and Treatment in a Cox proportional hazards model to assess survival. Information here: https://cran.r-project.org/web/packages/coxme/vignettes/coxme.pdf. “Surv” creates a survival object to combine the days column (NumDaysAlive) and the reproductive censor column (ReproductionCensor) to be used as a response variable in a model formula. The model follows the same syntax as linear models (lm, lmer, etc). Fixed effects: “var1 * var2” will give you the interaction term and the individual variables: so “var1 + var2 + var1:var2”. To add random effects, type “+ (1|random effect)”.
##Histograms
#Number of Days Alive - All data
hist(mydata$NumDaysAlive,col='steelblue',main='Original')
#Number of Days Alive - By Habitat
ggplot(mydata,aes(x=NumDaysAlive))+geom_histogram()+facet_wrap(vars(Habitat)) #Here we see that both Habitats have the same bi-modal distribution.
#Number of Days Alive - By Treatment
ggplot(mydata,aes(x=NumDaysAlive))+geom_histogram()+facet_wrap(vars(Treatment)) #Here we see that 4 of the treatments have the same bi-modal distribution and Biomass Removal is right skewed.
##Cox Model Start by making a simple model with no random effects. This will be compared to the full model with random effects.
cox_simplemodel1<-coxph(Surv(NumDaysAlive,
CensorReproduction)~
Habitat*Treatment,
data=mydata)
summary(cox_simplemodel1)
Call:
coxph(formula = Surv(NumDaysAlive, CensorReproduction) ~ Habitat *
Treatment, data = mydata)
n= 800, number of events= 253
coef exp(coef) se(coef) z Pr(>|z|)
HabitatVegetated -0.5980 0.5499 0.3308 -1.808 0.0706 .
TreatmentBiomass Removal 1.2401 3.4559 0.2771 4.475 7.65e-06 ***
TreatmentHemizonia 0.1563 1.1692 0.3364 0.465 0.6422
TreatmentRaking -0.8104 0.4447 0.4404 -1.840 0.0657 .
TreatmentRaking & Clipping 0.3718 1.4504 0.3304 1.125 0.2604
HabitatVegetated:TreatmentBiomass Removal 0.5848 1.7946 0.3772 1.550 0.1210
HabitatVegetated:TreatmentHemizonia 0.2799 1.3229 0.4740 0.590 0.5549
HabitatVegetated:TreatmentRaking 1.1347 3.1102 0.5861 1.936 0.0529 .
HabitatVegetated:TreatmentRaking & Clipping 0.8470 2.3327 0.4581 1.849 0.0645 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
exp(coef) exp(-coef) lower .95 upper .95
HabitatVegetated 0.5499 1.8185 0.2875 1.052
TreatmentBiomass Removal 3.4559 0.2894 2.0076 5.949
TreatmentHemizonia 1.1692 0.8553 0.6047 2.261
TreatmentRaking 0.4447 2.2489 0.1876 1.054
TreatmentRaking & Clipping 1.4504 0.6895 0.7590 2.771
HabitatVegetated:TreatmentBiomass Removal 1.7946 0.5572 0.8569 3.758
HabitatVegetated:TreatmentHemizonia 1.3229 0.7559 0.5224 3.350
HabitatVegetated:TreatmentRaking 3.1102 0.3215 0.9861 9.809
HabitatVegetated:TreatmentRaking & Clipping 2.3327 0.4287 0.9504 5.725
Concordance= 0.604 (se = 0.024 )
Likelihood ratio test= 90.68 on 9 df, p=1e-15
Wald test = 77.28 on 9 df, p=6e-13
Score (logrank) test = 87.13 on 9 df, p=6e-15
print(cox_simplemodel1)
Call:
coxph(formula = Surv(NumDaysAlive, CensorReproduction) ~ Habitat *
Treatment, data = mydata)
coef exp(coef) se(coef) z p
HabitatVegetated -0.5980 0.5499 0.3308 -1.808 0.0706
TreatmentBiomass Removal 1.2401 3.4559 0.2771 4.475 7.65e-06
TreatmentHemizonia 0.1563 1.1692 0.3364 0.465 0.6422
TreatmentRaking -0.8104 0.4447 0.4404 -1.840 0.0657
TreatmentRaking & Clipping 0.3718 1.4504 0.3304 1.125 0.2604
HabitatVegetated:TreatmentBiomass Removal 0.5848 1.7946 0.3772 1.550 0.1210
HabitatVegetated:TreatmentHemizonia 0.2799 1.3229 0.4740 0.590 0.5549
HabitatVegetated:TreatmentRaking 1.1347 3.1102 0.5861 1.936 0.0529
HabitatVegetated:TreatmentRaking & Clipping 0.8470 2.3327 0.4581 1.849 0.0645
Likelihood ratio test=90.68 on 9 df, p=1.187e-15
n= 800, number of events= 253
predict(cox_simplemodel1)
[1] 1.2400817 1.2400817 1.2400817 1.2400817 1.2400817 1.2400817 1.2400817 1.2400817
[9] 1.2400817 1.2400817 1.2400817 1.2400817 1.2400817 1.2400817 1.2400817 1.2400817
[17] 1.2400817 1.2400817 1.2400817 1.2400817 1.2400817 1.2400817 1.2400817 1.2400817
[25] 1.2400817 1.2400817 1.2400817 1.2400817 1.2400817 1.2400817 1.2400817 1.2400817
[33] 1.2400817 1.2400817 1.2400817 1.2400817 1.2400817 1.2400817 1.2400817 1.2400817
[41] 1.2400817 1.2400817 1.2400817 1.2400817 1.2400817 1.2400817 1.2400817 1.2400817
[49] 1.2400817 1.2400817 1.2400817 1.2400817 1.2400817 1.2400817 1.2400817 1.2400817
[57] 1.2400817 1.2400817 1.2400817 1.2400817 1.2400817 1.2400817 1.2400817 1.2400817
[65] 1.2400817 1.2400817 1.2400817 1.2400817 1.2400817 1.2400817 1.2400817 1.2400817
[73] 1.2400817 1.2400817 1.2400817 1.2400817 1.2400817 1.2400817 1.2400817 1.2400817
[81] 1.2268298 1.2268298 1.2268298 1.2268298 1.2268298 1.2268298 1.2268298 1.2268298
[89] 1.2268298 1.2268298 1.2268298 1.2268298 1.2268298 1.2268298 1.2268298 1.2268298
[97] 1.2268298 1.2268298 1.2268298 1.2268298 1.2268298 1.2268298 1.2268298 1.2268298
[105] 1.2268298 1.2268298 1.2268298 1.2268298 1.2268298 1.2268298 1.2268298 1.2268298
[113] 1.2268298 1.2268298 1.2268298 1.2268298 1.2268298 1.2268298 1.2268298 1.2268298
[121] 1.2268298 1.2268298 1.2268298 1.2268298 1.2268298 1.2268298 1.2268298 1.2268298
[129] 1.2268298 1.2268298 1.2268298 1.2268298 1.2268298 1.2268298 1.2268298 1.2268298
[137] 1.2268298 1.2268298 1.2268298 1.2268298 1.2268298 1.2268298 1.2268298 1.2268298
[145] 1.2268298 1.2268298 1.2268298 1.2268298 1.2268298 1.2268298 1.2268298 1.2268298
[153] 1.2268298 1.2268298 1.2268298 1.2268298 1.2268298 1.2268298 1.2268298 1.2268298
[161] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[169] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[177] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[185] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[193] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[201] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[209] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[217] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[225] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[233] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[241] -0.5980207 -0.5980207 -0.5980207 -0.5980207 -0.5980207 -0.5980207 -0.5980207 -0.5980207
[249] -0.5980207 -0.5980207 -0.5980207 -0.5980207 -0.5980207 -0.5980207 -0.5980207 -0.5980207
[257] -0.5980207 -0.5980207 -0.5980207 -0.5980207 -0.5980207 -0.5980207 -0.5980207 -0.5980207
[265] -0.5980207 -0.5980207 -0.5980207 -0.5980207 -0.5980207 -0.5980207 -0.5980207 -0.5980207
[273] -0.5980207 -0.5980207 -0.5980207 -0.5980207 -0.5980207 -0.5980207 -0.5980207 -0.5980207
[281] -0.5980207 -0.5980207 -0.5980207 -0.5980207 -0.5980207 -0.5980207 -0.5980207 -0.5980207
[289] -0.5980207 -0.5980207 -0.5980207 -0.5980207 -0.5980207 -0.5980207 -0.5980207 -0.5980207
[297] -0.5980207 -0.5980207 -0.5980207 -0.5980207 -0.5980207 -0.5980207 -0.5980207 -0.5980207
[305] -0.5980207 -0.5980207 -0.5980207 -0.5980207 -0.5980207 -0.5980207 -0.5980207 -0.5980207
[313] -0.5980207 -0.5980207 -0.5980207 -0.5980207 -0.5980207 -0.5980207 -0.5980207 -0.5980207
[321] 0.1563201 0.1563201 0.1563201 0.1563201 0.1563201 0.1563201 0.1563201 0.1563201
[329] 0.1563201 0.1563201 0.1563201 0.1563201 0.1563201 0.1563201 0.1563201 0.1563201
[337] 0.1563201 0.1563201 0.1563201 0.1563201 0.1563201 0.1563201 0.1563201 0.1563201
[345] 0.1563201 0.1563201 0.1563201 0.1563201 0.1563201 0.1563201 0.1563201 0.1563201
[353] 0.1563201 0.1563201 0.1563201 0.1563201 0.1563201 0.1563201 0.1563201 0.1563201
[361] 0.1563201 0.1563201 0.1563201 0.1563201 0.1563201 0.1563201 0.1563201 0.1563201
[369] 0.1563201 0.1563201 0.1563201 0.1563201 0.1563201 0.1563201 0.1563201 0.1563201
[377] 0.1563201 0.1563201 0.1563201 0.1563201 0.1563201 0.1563201 0.1563201 0.1563201
[385] 0.1563201 0.1563201 0.1563201 0.1563201 0.1563201 0.1563201 0.1563201 0.1563201
[393] 0.1563201 0.1563201 0.1563201 0.1563201 0.1563201 0.1563201 0.1563201 0.1563201
[401] -0.1618398 -0.1618398 -0.1618398 -0.1618398 -0.1618398 -0.1618398 -0.1618398 -0.1618398
[409] -0.1618398 -0.1618398 -0.1618398 -0.1618398 -0.1618398 -0.1618398 -0.1618398 -0.1618398
[417] -0.1618398 -0.1618398 -0.1618398 -0.1618398 -0.1618398 -0.1618398 -0.1618398 -0.1618398
[425] -0.1618398 -0.1618398 -0.1618398 -0.1618398 -0.1618398 -0.1618398 -0.1618398 -0.1618398
[433] -0.1618398 -0.1618398 -0.1618398 -0.1618398 -0.1618398 -0.1618398 -0.1618398 -0.1618398
[441] -0.1618398 -0.1618398 -0.1618398 -0.1618398 -0.1618398 -0.1618398 -0.1618398 -0.1618398
[449] -0.1618398 -0.1618398 -0.1618398 -0.1618398 -0.1618398 -0.1618398 -0.1618398 -0.1618398
[457] -0.1618398 -0.1618398 -0.1618398 -0.1618398 -0.1618398 -0.1618398 -0.1618398 -0.1618398
[465] -0.1618398 -0.1618398 -0.1618398 -0.1618398 -0.1618398 -0.1618398 -0.1618398 -0.1618398
[473] -0.1618398 -0.1618398 -0.1618398 -0.1618398 -0.1618398 -0.1618398 -0.1618398 -0.1618398
[481] -0.8104351 -0.8104351 -0.8104351 -0.8104351 -0.8104351 -0.8104351 -0.8104351 -0.8104351
[489] -0.8104351 -0.8104351 -0.8104351 -0.8104351 -0.8104351 -0.8104351 -0.8104351 -0.8104351
[497] -0.8104351 -0.8104351 -0.8104351 -0.8104351 -0.8104351 -0.8104351 -0.8104351 -0.8104351
[505] -0.8104351 -0.8104351 -0.8104351 -0.8104351 -0.8104351 -0.8104351 -0.8104351 -0.8104351
[513] -0.8104351 -0.8104351 -0.8104351 -0.8104351 -0.8104351 -0.8104351 -0.8104351 -0.8104351
[521] -0.8104351 -0.8104351 -0.8104351 -0.8104351 -0.8104351 -0.8104351 -0.8104351 -0.8104351
[529] -0.8104351 -0.8104351 -0.8104351 -0.8104351 -0.8104351 -0.8104351 -0.8104351 -0.8104351
[537] -0.8104351 -0.8104351 -0.8104351 -0.8104351 -0.8104351 -0.8104351 -0.8104351 -0.8104351
[545] -0.8104351 -0.8104351 -0.8104351 -0.8104351 -0.8104351 -0.8104351 -0.8104351 -0.8104351
[553] -0.8104351 -0.8104351 -0.8104351 -0.8104351 -0.8104351 -0.8104351 -0.8104351 -0.8104351
[561] -0.2737686 -0.2737686 -0.2737686 -0.2737686 -0.2737686 -0.2737686 -0.2737686 -0.2737686
[569] -0.2737686 -0.2737686 -0.2737686 -0.2737686 -0.2737686 -0.2737686 -0.2737686 -0.2737686
[577] -0.2737686 -0.2737686 -0.2737686 -0.2737686 -0.2737686 -0.2737686 -0.2737686 -0.2737686
[585] -0.2737686 -0.2737686 -0.2737686 -0.2737686 -0.2737686 -0.2737686 -0.2737686 -0.2737686
[593] -0.2737686 -0.2737686 -0.2737686 -0.2737686 -0.2737686 -0.2737686 -0.2737686 -0.2737686
[601] -0.2737686 -0.2737686 -0.2737686 -0.2737686 -0.2737686 -0.2737686 -0.2737686 -0.2737686
[609] -0.2737686 -0.2737686 -0.2737686 -0.2737686 -0.2737686 -0.2737686 -0.2737686 -0.2737686
[617] -0.2737686 -0.2737686 -0.2737686 -0.2737686 -0.2737686 -0.2737686 -0.2737686 -0.2737686
[625] -0.2737686 -0.2737686 -0.2737686 -0.2737686 -0.2737686 -0.2737686 -0.2737686 -0.2737686
[633] -0.2737686 -0.2737686 -0.2737686 -0.2737686 -0.2737686 -0.2737686 -0.2737686 -0.2737686
[641] 0.3718368 0.3718368 0.3718368 0.3718368 0.3718368 0.3718368 0.3718368 0.3718368
[649] 0.3718368 0.3718368 0.3718368 0.3718368 0.3718368 0.3718368 0.3718368 0.3718368
[657] 0.3718368 0.3718368 0.3718368 0.3718368 0.3718368 0.3718368 0.3718368 0.3718368
[665] 0.3718368 0.3718368 0.3718368 0.3718368 0.3718368 0.3718368 0.3718368 0.3718368
[673] 0.3718368 0.3718368 0.3718368 0.3718368 0.3718368 0.3718368 0.3718368 0.3718368
[681] 0.3718368 0.3718368 0.3718368 0.3718368 0.3718368 0.3718368 0.3718368 0.3718368
[689] 0.3718368 0.3718368 0.3718368 0.3718368 0.3718368 0.3718368 0.3718368 0.3718368
[697] 0.3718368 0.3718368 0.3718368 0.3718368 0.3718368 0.3718368 0.3718368 0.3718368
[705] 0.3718368 0.3718368 0.3718368 0.3718368 0.3718368 0.3718368 0.3718368 0.3718368
[713] 0.3718368 0.3718368 0.3718368 0.3718368 0.3718368 0.3718368 0.3718368 0.3718368
[721] 0.6208270 0.6208270 0.6208270 0.6208270 0.6208270 0.6208270 0.6208270 0.6208270
[729] 0.6208270 0.6208270 0.6208270 0.6208270 0.6208270 0.6208270 0.6208270 0.6208270
[737] 0.6208270 0.6208270 0.6208270 0.6208270 0.6208270 0.6208270 0.6208270 0.6208270
[745] 0.6208270 0.6208270 0.6208270 0.6208270 0.6208270 0.6208270 0.6208270 0.6208270
[753] 0.6208270 0.6208270 0.6208270 0.6208270 0.6208270 0.6208270 0.6208270 0.6208270
[761] 0.6208270 0.6208270 0.6208270 0.6208270 0.6208270 0.6208270 0.6208270 0.6208270
[769] 0.6208270 0.6208270 0.6208270 0.6208270 0.6208270 0.6208270 0.6208270 0.6208270
[777] 0.6208270 0.6208270 0.6208270 0.6208270 0.6208270 0.6208270 0.6208270 0.6208270
[785] 0.6208270 0.6208270 0.6208270 0.6208270 0.6208270 0.6208270 0.6208270 0.6208270
[793] 0.6208270 0.6208270 0.6208270 0.6208270 0.6208270 0.6208270 0.6208270 0.6208270
hist(predict(cox_simplemodel1))
Now make a full model using random effects
cox_fullmodel1<-coxme(Surv(NumDaysAlive,CensorReproduction)~
Habitat*Treatment+
(1|Site)+
(1|Block),
data=mydata)
summary(cox_fullmodel1)
Cox mixed-effects model fit by maximum likelihood
Data: mydata
events, n = 253, 800
Iterations= 19 119
NULL Integrated Fitted
Log-likelihood -1204.23 -1117.392 -1093.286
Chisq df p AIC BIC
Integrated loglik 173.68 11.00 0 151.68 112.81
Penalized loglik 221.89 22.42 0 177.04 97.81
Model: Surv(NumDaysAlive, CensorReproduction) ~ Habitat * Treatment + (1 | Site) + (1 | Block)
Fixed coefficients
coef exp(coef) se(coef) z p
HabitatVegetated -0.37919954 0.6844090 0.3444499 -1.10 2.7e-01
TreatmentBiomass Removal 1.38647395 4.0007184 0.2911902 4.76 1.9e-06
TreatmentHemizonia -0.09751812 0.9070859 0.3540413 -0.28 7.8e-01
TreatmentRaking -0.64329505 0.5255578 0.4572637 -1.41 1.6e-01
TreatmentRaking & Clipping 0.59622662 1.8152562 0.3498383 1.70 8.8e-02
HabitatVegetated:TreatmentBiomass Removal 0.45521372 1.5765103 0.3874747 1.17 2.4e-01
HabitatVegetated:TreatmentHemizonia 0.25268759 1.2874810 0.4903683 0.52 6.1e-01
HabitatVegetated:TreatmentRaking 0.80143838 2.2287444 0.6000419 1.34 1.8e-01
HabitatVegetated:TreatmentRaking & Clipping 0.68420207 1.9821896 0.4766792 1.44 1.5e-01
Random effects
Group Variable Std Dev Variance
Site Intercept 0.29585965 0.08753293
Block Intercept 0.83465076 0.69664189
print(cox_fullmodel1)
Cox mixed-effects model fit by maximum likelihood
Data: mydata
events, n = 253, 800
Iterations= 19 119
NULL Integrated Fitted
Log-likelihood -1204.23 -1117.392 -1093.286
Chisq df p AIC BIC
Integrated loglik 173.68 11.00 0 151.68 112.81
Penalized loglik 221.89 22.42 0 177.04 97.81
Model: Surv(NumDaysAlive, CensorReproduction) ~ Habitat * Treatment + (1 | Site) + (1 | Block)
Fixed coefficients
coef exp(coef) se(coef) z p
HabitatVegetated -0.37919954 0.6844090 0.3444499 -1.10 2.7e-01
TreatmentBiomass Removal 1.38647395 4.0007184 0.2911902 4.76 1.9e-06
TreatmentHemizonia -0.09751812 0.9070859 0.3540413 -0.28 7.8e-01
TreatmentRaking -0.64329505 0.5255578 0.4572637 -1.41 1.6e-01
TreatmentRaking & Clipping 0.59622662 1.8152562 0.3498383 1.70 8.8e-02
HabitatVegetated:TreatmentBiomass Removal 0.45521372 1.5765103 0.3874747 1.17 2.4e-01
HabitatVegetated:TreatmentHemizonia 0.25268759 1.2874810 0.4903683 0.52 6.1e-01
HabitatVegetated:TreatmentRaking 0.80143838 2.2287444 0.6000419 1.34 1.8e-01
HabitatVegetated:TreatmentRaking & Clipping 0.68420207 1.9821896 0.4766792 1.44 1.5e-01
Random effects
Group Variable Std Dev Variance
Site Intercept 0.29585965 0.08753293
Block Intercept 0.83465076 0.69664189
predict(cox_fullmodel1)
[1] 3.412069248 1.723904182 1.142354944 1.411114662 2.826926974 1.710423900
[7] 1.521797855 1.543482099 1.020912356 0.832419259 3.057516702 1.369351636
[13] 0.787802397 1.056562116 2.472374428 1.355871354 1.167245309 1.188929553
[19] 0.666359810 0.477866713 3.027300084 1.339135018 0.757585779 1.026345498
[25] 2.442157809 1.325654736 1.137028691 1.158712935 0.636143192 0.447650095
[31] 3.111607266 1.423442200 0.841892961 1.110652680 2.526464992 1.409961918
[37] 1.221335873 1.243020117 0.720450374 0.531957277 2.783497956 1.095332890
[43] 0.513783652 0.782543370 2.198355682 1.081852608 0.893226563 0.914910807
[49] 0.392341064 0.203847967 3.491439676 1.803274610 1.221725371 1.490485090
[55] 2.906297402 1.789794328 1.601168283 1.622852527 1.100282784 0.911789687
[61] 2.874174708 1.186009641 0.604460403 0.873220121 2.289032433 1.172529359
[67] 0.983903315 1.005587559 0.483017815 0.294524719 2.914415579 1.226250512
[73] 0.644701274 0.913460992 2.329273304 1.212770231 1.024144186 1.045828430
[79] 0.523258686 0.334765590 3.488083434 1.799918367 1.218369129 1.487128847
[85] 2.902941159 1.786438086 1.597812041 1.619496285 1.096926541 0.908433445
[91] 3.133530888 1.445365821 0.863816583 1.132576301 2.548388613 1.431885540
[97] 1.243259495 1.264943739 0.742373995 0.553880899 3.103314270 1.415149203
[103] 0.833599965 1.102359683 2.518171995 1.401668921 1.213042877 1.234727121
[109] 0.712157377 0.523664281 3.187621452 1.499456385 0.917907147 1.186666865
[115] 2.602479177 1.485976104 1.297350059 1.319034303 0.796464559 0.607971463
[121] 2.859512142 1.171347075 0.589797837 0.858557555 2.274369867 1.157866794
[127] 0.969240749 0.990924993 0.468355249 0.279862153 3.567453862 1.879288795
[133] 1.297739557 1.566499275 2.982311587 1.865808513 1.677182469 1.698866713
[139] 1.176296969 0.987803873 2.950188893 1.262023827 0.680474588 0.949234307
[145] 2.365046619 1.248543545 1.059917500 1.081601744 0.559032001 0.370538904
[151] 2.990429765 1.302264698 0.720715460 0.989475178 2.405287490 1.288784416
[157] 1.100158372 1.121842616 0.599272872 0.410779776 2.025595296 0.337430230
[163] -0.244119009 0.024640710 1.440453022 0.323949948 0.135323903 0.157008147
[169] -0.365561596 -0.554054693 1.671042750 -0.017122316 -0.598671555 -0.329911836
[175] 1.085900476 -0.030602598 -0.219228643 -0.197544399 -0.720114142 -0.908607239
[181] 1.640826132 -0.047338935 -0.628888173 -0.360128455 1.055683857 -0.060819216
[187] -0.249445261 -0.227761017 -0.750330761 -0.938823857 1.725133314 0.036968248
[193] -0.544580991 -0.275821272 1.139991040 0.023487966 -0.165138079 -0.143453835
[199] -0.666023578 -0.854516675 1.397024004 -0.291141062 -0.872690301 -0.603930582
[205] 0.811881730 -0.304621344 -0.493247389 -0.471563145 -0.994132888 -1.182625985
[211] 2.104965724 0.416800657 -0.164748581 0.104011137 1.519823449 0.403320376
[217] 0.214694331 0.236378575 -0.286191169 -0.474684265 1.487700756 -0.200464311
[223] -0.782013549 -0.513253831 0.902558481 -0.213944593 -0.402570638 -0.380886393
[229] -0.903456137 -1.091949233 1.527941627 -0.160223440 -0.741772678 -0.473012960
[235] 0.942799352 -0.173703721 -0.362329766 -0.340645522 -0.863215266 -1.051708362
[241] 1.646395758 -0.041769308 -0.623318547 -0.354558828 1.061253484 -0.055249590
[247] -0.243875635 -0.222191391 -0.744761134 -0.933254231 1.291843212 -0.396321855
[253] -0.977871093 -0.709111374 0.706700937 -0.409802136 -0.598428181 -0.576743937
[259] -1.099313681 -1.287806777 1.261626594 -0.426538473 -1.008087711 -0.739327993
[265] 0.676484319 -0.440018754 -0.628644799 -0.606960555 -1.129530299 -1.318023395
[271] 1.345933776 -0.342231290 -0.923780529 -0.655020810 0.760791502 -0.355711572
[277] -0.544337617 -0.522653373 -1.045223116 -1.233716213 1.017824466 -0.670340600
[283] -1.251889839 -0.983130120 0.432682192 -0.683820882 -0.872446927 -0.850762683
[289] -1.373332426 -1.561825523 1.725766186 0.037601119 -0.543948119 -0.275188401
[295] 1.140623911 0.024120838 -0.164505207 -0.142820963 -0.665390707 -0.853883803
[301] 1.108501217 -0.579663849 -1.161213087 -0.892453369 0.523358943 -0.593144131
[307] -0.781770176 -0.760085932 -1.282655675 -1.471148772 1.148742089 -0.539422978
[313] -1.120972216 -0.852212498 0.563599814 -0.552903260 -0.741529304 -0.719845060
[319] -1.242414804 -1.430907900 1.928077178 0.239912111 -0.341637127 -0.072877409
[325] 1.342934903 0.226431829 0.037805785 0.059490029 -0.463079715 -0.651572811
[331] 1.573524632 -0.114640435 -0.696189673 -0.427429955 0.988382357 -0.128120717
[337] -0.316746761 -0.295062517 -0.817632261 -1.006125357 1.543308014 -0.144857053
[343] -0.726406291 -0.457646573 0.958165739 -0.158337335 -0.346963380 -0.325279135
[349] -0.847848879 -1.036341975 1.627615196 -0.060549871 -0.642099109 -0.373339391
[355] 1.042472921 -0.074030153 -0.262656197 -0.240971953 -0.763541697 -0.952034793
[361] 1.299505886 -0.388659181 -0.970208419 -0.701448701 0.714363611 -0.402139463
[367] -0.590765507 -0.569081263 -1.091651007 -1.280144103 2.007447606 0.319282539
[373] -0.262266699 0.006493019 1.422305331 0.305802257 0.117176212 0.138860457
[379] -0.383709287 -0.572202383 1.390182637 -0.297982430 -0.879531668 -0.610771950
[385] 0.805040362 -0.311462711 -0.500088756 -0.478404512 -1.000974256 -1.189467352
[391] 1.430423508 -0.257741558 -0.839290797 -0.570531078 0.845281234 -0.271221840
[397] -0.459847885 -0.438163641 -0.960733384 -1.149226481 1.801565225 0.113400159
[403] -0.468149079 -0.199389361 1.216422951 0.099919877 -0.088706168 -0.067021924
[409] -0.589591667 -0.778084764 1.447012679 -0.241152387 -0.822701626 -0.553941907
[415] 0.861870405 -0.254632669 -0.443258714 -0.421574470 -0.944144213 -1.132637310
[421] 1.416796061 -0.271369005 -0.852918244 -0.584158525 0.831653787 -0.284849287
[427] -0.473475332 -0.451791088 -0.974360831 -1.162853928 1.501103243 -0.187061823
[433] -0.768611061 -0.499851343 0.915960969 -0.200542105 -0.389168150 -0.367483906
[439] -0.890053649 -1.078546746 1.172993933 -0.515171133 -1.096720371 -0.827960653
[445] 0.587851659 -0.528651415 -0.717277460 -0.695593216 -1.218162959 -1.406656056
[451] 1.880935653 0.192770587 -0.388778652 -0.120018933 1.295793379 0.179290305
[457] -0.009335740 0.012348504 -0.510221239 -0.698714336 1.263670685 -0.424494382
[463] -1.006043620 -0.737283902 0.678528410 -0.437974664 -0.626600708 -0.604916464
[469] -1.127486208 -1.315979304 1.303911556 -0.384253511 -0.965802749 -0.697043031
[475] 0.718769281 -0.397733792 -0.586359837 -0.564675593 -1.087245337 -1.275738433
[481] 1.382300247 -0.305864820 -0.887414058 -0.618654340 0.797157972 -0.319345102
[487] -0.507971146 -0.486286902 -1.008856646 -1.197349742 1.027747700 -0.660417366
[493] -1.241966604 -0.973206886 0.442605426 -0.673897648 -0.862523693 -0.840839449
[499] -1.363409192 -1.551902288 0.997531082 -0.690633984 -1.272183223 -1.003423504
[505] 0.412388808 -0.704114266 -0.892740311 -0.871056067 -1.393625810 -1.582118907
[511] 1.081838265 -0.606326802 -1.187876040 -0.919116322 0.496695990 -0.619807084
[517] -0.808433128 -0.786748884 -1.309318628 -1.497811724 0.753728955 -0.934436112
[523] -1.515985350 -1.247225632 0.168586680 -0.947916394 -1.136542438 -1.114858194
[529] -1.637427938 -1.825921034 1.461670674 -0.226494392 -0.808043631 -0.539283912
[535] 0.876528400 -0.239974674 -0.428600719 -0.406916475 -0.929486218 -1.117979315
[541] 0.844405706 -0.843759361 -1.425308599 -1.156548881 0.259263431 -0.857239643
[547] -1.045865687 -1.024181443 -1.546751187 -1.735244283 0.884646577 -0.803518489
[553] -1.385067728 -1.116308009 0.299504303 -0.816998771 -1.005624816 -0.983940572
[559] -1.506510315 -1.695003412 1.804539087 0.116374020 -0.465175218 -0.196415500
[565] 1.219396812 0.102893738 -0.085732307 -0.064048062 -0.586617806 -0.775110902
[571] 1.449986540 -0.238178526 -0.819727765 -0.550968046 0.864844266 -0.251658808
[577] -0.440284853 -0.418600609 -0.941170352 -1.129663449 1.419769922 -0.268395144
[583] -0.849944383 -0.581184664 0.834627648 -0.281875426 -0.470501471 -0.448817227
[589] -0.971386970 -1.159880067 1.504077104 -0.184087962 -0.765637200 -0.496877482
[595] 0.918934830 -0.197568244 -0.386194289 -0.364510045 -0.887079788 -1.075572884
[601] 1.175967795 -0.512197272 -1.093746510 -0.824986792 0.590825520 -0.525677554
[607] -0.714303599 -0.692619355 -1.215189098 -1.403682194 1.883909514 0.195744448
[613] -0.385804791 -0.117045072 1.298767240 0.182264166 -0.006361879 0.015322365
[619] -0.507247378 -0.695740475 1.266644546 -0.421520521 -1.003069759 -0.734310041
[625] 0.681502271 -0.435000803 -0.623626847 -0.601942603 -1.124512347 -1.313005443
[631] 1.306885417 -0.381279650 -0.962828888 -0.694069170 0.721743142 -0.394759931
[637] -0.583385976 -0.561701732 -1.084271476 -1.272764572 2.621821912 0.933656845
[643] 0.352107607 0.620867325 2.036679637 0.920176564 0.731550519 0.753234763
[649] 0.230665019 0.042171923 2.267269366 0.579104299 -0.002444939 0.266314779
[655] 1.682127091 0.565624018 0.376997973 0.398682217 -0.123887527 -0.312380623
[661] 2.237052748 0.548887681 -0.032661557 0.236098161 1.651910473 0.535407399
[667] 0.346781355 0.368465599 -0.154104145 -0.342597241 2.321359930 0.633194863
[673] 0.051645625 0.320405343 1.736217655 0.619714582 0.431088537 0.452772781
[679] -0.069796963 -0.258290059 1.993250620 0.305085553 -0.276463685 -0.007703967
[685] 1.408108345 0.291605272 0.102979227 0.124663471 -0.397906273 -0.586399369
[691] 2.701192340 1.013027273 0.431478035 0.700237753 2.116050065 0.999546991
[697] 0.810920947 0.832605191 0.310035447 0.121542351 2.083927371 0.395762305
[703] -0.185786934 0.082972785 1.498785097 0.382282023 0.193655978 0.215340222
[709] -0.307229521 -0.495722618 2.124168243 0.436003176 -0.145546062 0.123213656
[715] 1.539025968 0.422522894 0.233896850 0.255581094 -0.266988650 -0.455481746
[721] 2.926824446 1.238659380 0.657110141 0.925869860 2.341682172 1.225179098
[727] 1.036553053 1.058237297 0.535667554 0.347174457 2.572271900 0.884106833
[733] 0.302557595 0.571317313 1.987129625 0.870626552 0.682000507 0.703684751
[739] 0.181115007 -0.007378089 2.542055282 0.853890215 0.272340977 0.541100695
[745] 1.956913007 0.840409934 0.651783889 0.673468133 0.150898389 -0.037594707
[751] 2.626362464 0.938197398 0.356648159 0.625407878 2.041220189 0.924717116
[757] 0.736091071 0.757775315 0.235205572 0.046712475 2.298253154 0.610088088
[763] 0.028538849 0.297298568 1.713110880 0.596607806 0.407981761 0.429666005
[769] -0.092903738 -0.281396835 3.006194874 1.318029807 0.736480569 1.005240287
[775] 2.421052599 1.304549526 1.115923481 1.137607725 0.615037981 0.426544885
[781] 2.388929905 0.700764839 0.119215600 0.387975319 1.803787631 0.687284557
[787] 0.498658512 0.520342756 -0.002226987 -0.190720084 2.429170777 0.741005710
[793] 0.159456472 0.428216190 1.844028502 0.727525428 0.538899384 0.560583628
[799] 0.038013884 -0.150479212
hist(predict(cox_fullmodel1))
Now we can compare the models to see which model is best
anova(cox_simplemodel1,cox_fullmodel1)
Analysis of Deviance Table
Cox model: response is Surv(NumDaysAlive, CensorReproduction)
Model 1: ~Habitat * Treatment
Model 2: ~Habitat * Treatment + (1 | Site) + (1 | Block)
loglik Chisq Df P(>|Chi|)
1 -1158.9
2 -1117.4 82.993 2 < 2.2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#See example: https://www.rdocumentation.org/packages/coxme/versions/2.2-16/topics/coxme
AIC(cox_simplemodel1,cox_fullmodel1) #But comparing the AIC scores is easiest. Keep the lower AIC score because that is considered the better model. Here it is the fullmodel1.
Now, let’s make a model with no interaction term
cox_fullmodel2<-coxme(Surv(NumDaysAlive,CensorReproduction)~
Habitat+
Treatment+
(1|Site)+
(1|Block),
data=mydata)
summary(cox_fullmodel2)
Cox mixed-effects model fit by maximum likelihood
Data: mydata
events, n = 253, 800
Iterations= 18 95
NULL Integrated Fitted
Log-likelihood -1204.23 -1118.851 -1094.381
Chisq df p AIC BIC
Integrated loglik 170.76 7.00 0 156.76 132.02
Penalized loglik 219.70 18.58 0 182.53 116.86
Model: Surv(NumDaysAlive, CensorReproduction) ~ Habitat + Treatment + (1 | Site) + (1 | Block)
Fixed coefficients
coef exp(coef) se(coef) z p
HabitatVegetated 0.04968034 1.0509351 0.1326079 0.37 7.1e-01
TreatmentBiomass Removal 1.60910167 4.9983191 0.2268731 7.09 1.3e-12
TreatmentHemizonia 0.01272746 1.0128088 0.2465085 0.05 9.6e-01
TreatmentRaking -0.20667151 0.8132868 0.3019860 -0.68 4.9e-01
TreatmentRaking & Clipping 0.95534162 2.5995585 0.2426422 3.94 8.2e-05
Random effects
Group Variable Std Dev Variance
Site Intercept 0.30795467 0.09483608
Block Intercept 0.84503192 0.71407895
Let’s compare the the first two models to test for the significance of the term that is removed (using LR)
anova(cox_fullmodel1,cox_fullmodel2) #Not significant
Analysis of Deviance Table
Cox model: response is Surv(NumDaysAlive, CensorReproduction)
Model 1: ~Habitat * Treatment + (1 | Site) + (1 | Block)
Model 2: ~Habitat + Treatment + (1 | Site) + (1 | Block)
loglik Chisq Df P(>|Chi|)
1 -1117.4
2 -1118.8 2.919 4 0.5715
AIC(cox_fullmodel1,cox_fullmodel2) #Interaction term is not significant and the second model has a lower AIC score. So we can drop the interaction term and keep fullmodel2.
So, now let’s add in population nested under site as a random effect
cox_fullmodel3<-coxme(Surv(NumDaysAlive,CensorReproduction)~
Habitat+
Treatment+
(1|Site/Population)+
(1|Block),
data=mydata)
summary(cox_fullmodel3)
Cox mixed-effects model fit by maximum likelihood
Data: mydata
events, n = 253, 800
Iterations= 20 105
NULL Integrated Fitted
Log-likelihood -1204.23 -1116.667 -1086.004
Chisq df p AIC BIC
Integrated loglik 175.13 8.00 0 159.13 130.86
Penalized loglik 236.45 22.68 0 191.10 110.96
Model: Surv(NumDaysAlive, CensorReproduction) ~ Habitat + Treatment + (1 | Site/Population) + (1 | Block)
Fixed coefficients
coef exp(coef) se(coef) z p
HabitatVegetated 0.03971413 1.0405133 0.2256647 0.18 8.6e-01
TreatmentBiomass Removal 1.58216404 4.8654735 0.2279166 6.94 3.9e-12
TreatmentHemizonia -0.03665974 0.9640041 0.2497694 -0.15 8.8e-01
TreatmentRaking -0.27759259 0.7576054 0.3053999 -0.91 3.6e-01
TreatmentRaking & Clipping 0.94413698 2.5705939 0.2466166 3.83 1.3e-04
Random effects
Group Variable Std Dev Variance
Site/Population (Intercept) 0.3669601665 0.1346597638
Site (Intercept) 0.0200515483 0.0004020646
Block Intercept 0.8757403665 0.7669211894
Now we can compare the models to see which model is best
anova(cox_fullmodel2,cox_fullmodel3) #Significant
Analysis of Deviance Table
Cox model: response is Surv(NumDaysAlive, CensorReproduction)
Model 1: ~Habitat + Treatment + (1 | Site) + (1 | Block)
Model 2: ~Habitat + Treatment + (1 | Site/Population) + (1 | Block)
loglik Chisq Df P(>|Chi|)
1 -1118.8
2 -1116.7 4.3685 1 0.03661 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
AIC(cox_fullmodel2,cox_fullmodel3) #Looks like fullmodel3 is the better model because of the lower AIC score
###Best Model
cox_fullmodel3<-coxme(Surv(NumDaysAlive,CensorReproduction)~
Habitat+
Treatment+
(1|Site/Population)+
(1|Block),
data=mydata)
summary(cox_fullmodel3)
Cox mixed-effects model fit by maximum likelihood
Data: mydata
events, n = 253, 800
Iterations= 20 105
NULL Integrated Fitted
Log-likelihood -1204.23 -1116.667 -1086.004
Chisq df p AIC BIC
Integrated loglik 175.13 8.00 0 159.13 130.86
Penalized loglik 236.45 22.68 0 191.10 110.96
Model: Surv(NumDaysAlive, CensorReproduction) ~ Habitat + Treatment + (1 | Site/Population) + (1 | Block)
Fixed coefficients
coef exp(coef) se(coef) z p
HabitatVegetated 0.03971413 1.0405133 0.2256647 0.18 8.6e-01
TreatmentBiomass Removal 1.58216404 4.8654735 0.2279166 6.94 3.9e-12
TreatmentHemizonia -0.03665974 0.9640041 0.2497694 -0.15 8.8e-01
TreatmentRaking -0.27759259 0.7576054 0.3053999 -0.91 3.6e-01
TreatmentRaking & Clipping 0.94413698 2.5705939 0.2466166 3.83 1.3e-04
Random effects
Group Variable Std Dev Variance
Site/Population (Intercept) 0.3669601665 0.1346597638
Site (Intercept) 0.0200515483 0.0004020646
Block Intercept 0.8757403665 0.7669211894
###Risk Assessment These values are found in the model summary, but if you want to pull them out, here is how you interpret them. 1 = no effect, <1 = decreased risk of death, >1 = increased risk of death.
exp(coef(cox_fullmodel3)) #This should be interpreted that Biomass Removal is almost 5% more likely to survive to reproduction compared to Control, and Raking + Clipping is about 2.5% more likely to survive to reproduction compared to Control.
HabitatVegetated TreatmentBiomass Removal TreatmentHemizonia
1.0405133 4.8654735 0.9640041
TreatmentRaking TreatmentRaking & Clipping
0.7576054 2.5705939
exp(ranef(cox_fullmodel3)$Block)
A B C D E F G H I
5.9382701 1.0302514 0.5553219 0.7170038 3.1240541 1.0901227 0.7725668 0.8326539 0.4957931
J
0.3779511
exp(ranef(cox_fullmodel3)$Site) # Pretty even among all the sites
BAY CHE GUA LEX OAP PAR PEN SSJ
1.0018910 0.9997332 0.9997487 1.0001393 0.9987233 1.0021430 0.9986960 0.9989318
Anova(cox_fullmodel3)
Analysis of Deviance Table (Type II tests)
Response: Surv(NumDaysAlive, CensorReproduction)
Df Chisq Pr(>Chisq)
Habitat 1 0.031 0.8603
Treatment 4 84.060 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Looks like roadside and offroad plants are the same, so let’s combine them together in a model (aka, removing the Habitat term)
cox_fullmodel4<-coxme(Surv(NumDaysAlive,CensorReproduction)~
Treatment+
(1|Site/Population)+
(1|Block),
data=mydata)
summary(cox_fullmodel4)
Cox mixed-effects model fit by maximum likelihood
Data: mydata
events, n = 253, 800
Iterations= 22 115
NULL Integrated Fitted
Log-likelihood -1204.23 -1116.683 -1086.019
Chisq df p AIC BIC
Integrated loglik 175.09 7.00 0 161.09 136.36
Penalized loglik 236.42 22.35 0 191.73 112.77
Model: Surv(NumDaysAlive, CensorReproduction) ~ Treatment + (1 | Site/Population) + (1 | Block)
Fixed coefficients
coef exp(coef) se(coef) z p
TreatmentBiomass Removal 1.5817104 4.8632667 0.2278906 6.94 3.9e-12
TreatmentHemizonia -0.0358444 0.9647904 0.2497418 -0.14 8.9e-01
TreatmentRaking -0.2777768 0.7574659 0.3054173 -0.91 3.6e-01
TreatmentRaking & Clipping 0.9445967 2.5717760 0.2466066 3.83 1.3e-04
Random effects
Group Variable Std Dev Variance
Site/Population (Intercept) 0.3670944471 0.1347583331
Site (Intercept) 0.0200501745 0.0004020095
Block Intercept 0.8752789365 0.7661132166
summary(cox_fullmodel3)
Cox mixed-effects model fit by maximum likelihood
Data: mydata
events, n = 253, 800
Iterations= 20 105
NULL Integrated Fitted
Log-likelihood -1204.23 -1116.667 -1086.004
Chisq df p AIC BIC
Integrated loglik 175.13 8.00 0 159.13 130.86
Penalized loglik 236.45 22.68 0 191.10 110.96
Model: Surv(NumDaysAlive, CensorReproduction) ~ Habitat + Treatment + (1 | Site/Population) + (1 | Block)
Fixed coefficients
coef exp(coef) se(coef) z p
HabitatVegetated 0.03971413 1.0405133 0.2256647 0.18 8.6e-01
TreatmentBiomass Removal 1.58216404 4.8654735 0.2279166 6.94 3.9e-12
TreatmentHemizonia -0.03665974 0.9640041 0.2497694 -0.15 8.8e-01
TreatmentRaking -0.27759259 0.7576054 0.3053999 -0.91 3.6e-01
TreatmentRaking & Clipping 0.94413698 2.5705939 0.2466166 3.83 1.3e-04
Random effects
Group Variable Std Dev Variance
Site/Population (Intercept) 0.3669601665 0.1346597638
Site (Intercept) 0.0200515483 0.0004020646
Block Intercept 0.8757403665 0.7669211894
anova(cox_fullmodel3,cox_fullmodel4) #Significant
Analysis of Deviance Table
Cox model: response is Surv(NumDaysAlive, CensorReproduction)
Model 1: ~Habitat + Treatment + (1 | Site/Population) + (1 | Block)
Model 2: ~Treatment + (1 | Site/Population) + (1 | Block)
loglik Chisq Df P(>|Chi|)
1 -1116.7
2 -1116.7 0.0319 1 0.8582
AIC(cox_fullmodel3,cox_fullmodel4) #Looks like fullmodel4 is the better model because the AIC score is within 2 points of each other, therefore the models are assessed the same and you should take the simpler model. But this is for another manuscript, probably.
##ggplot ###Interaction Plot
pd<-position_dodge(0)
surv_mydata1<-mydata%>%
filter(Treatment=="Biomass Removal"|Treatment=="Grassland")%>%
group_by(Treatment,Habitat)%>%
summarise(Mean=mean(PropBudSite),
SD=sd(PropBudSite),
N=length(PropBudSite))%>%
mutate(SE=SD/sqrt(N))
`summarise()` has grouped output by 'Treatment'. You can override using the `.groups` argument.
field_int_surv_all_gg<-ggplot(data=surv_mydata1,
aes(x=factor(Treatment,levels=c("Biomass Removal","Grassland")),
y=Mean,
shape=Habitat,
group=Habitat,
color=Habitat))+
ylab("Proportion Survival to Bud")+
coord_cartesian(ylim=c(0,1))+
xlab("")+
theme_classic(base_size=16)+
theme(panel.border=element_blank(),
panel.grid.major=element_blank(),
panel.grid.minor=element_blank(),
axis.line=element_line(colour="black"),
legend.position = c(0.8, 0.8))+
scale_y_continuous(name="Proportion Survival to Bud",
limits=c(0,1),
breaks=c(0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0))+
geom_line(size=1,
position=pd)+
geom_point(size=5,
position=pd)+
scale_color_manual(values=c("gray85","forestgreen"))+
geom_errorbar(width=0.15,
position=pd,
aes(ymin=(Mean-SE),
ymax=(Mean+SE)))+
labs(y="Proportion Survival to Bud",
fill="Source Habitat",
color="Source Habitat",
shape="Source Habitat")
field_int_surv_all_gg
ggsave(plot=field_int_surv_all_gg,file="/Users/Miranda/Documents/Education/UC Santa Cruz/Dittrichia/Dittrichia_Analysis/figures/field_int_surv_all_gg.png",width=25,height=15,units="cm",dpi=800)
##autoplot - Survival Curves Plot by habitat test
pd<-position_dodge(0)
cox_model_data1<-mydata%>%
filter(Treatment=="Biomass Removal"|Treatment=="Grassland")
cox_model_graph1<-survfit(Surv(NumDaysAlive,CensorReproduction)~
Habitat,
data=cox_model_data1)
surv_curve_habitat_treat<-ggplot(data=cox_model_graph1)+
aes(x=time,
y=surv,
shape=strata,
group=strata,
color=strata)+
theme_classic(base_size=16)+
scale_y_continuous(name="Survival Probability",
limits=c(0,1),
breaks=seq(0.0,1.0,0.1))+
scale_x_continuous(name="Survival Time (Days)",
limits=c(0,275),
breaks=seq(0,275,20))+
geom_step(size=1,
position=pd)+
geom_point(size=5,
position=pd)+
scale_color_manual(values=c("gray85","forestgreen"))+
labs(x="Survival Time (Days)",
y="Survival Probability")
surv_curve_habitat_treat
##autoplot - Survival Curves Plot by habitat
cox_model_graph1<-survfit(Surv(NumDaysAlive,CensorReproduction)~
Habitat,
data=mydata)
autoplot(cox_model_graph1)+
labs(x="\n Survival Time (Days)",
y="Survival Probabilities\n",
title="Survival Times Of \n Roadside and Vegetated Populations\n")+
theme_bw(base_size=16)+
theme(plot.title=element_text(hjust=0.5),
axis.title.x=element_text(face="bold",
color="Black",
size=12),
axis.title.y=element_text(face="bold",
color="Black",
size=12),
legend.title=element_text(face="bold",
size=10))
Plot by all treatment
cox_model_graph2<-survfit(Surv(NumDaysAlive,CensorReproduction)~
Treatment,
data=mydata)
autoplot(cox_model_graph2)+
labs(x="\n Survival Time (Days)",
y="Survival Probabilities\n",
title="Survival Times Of \n Roadside and Vegetated Populations\n")+
theme(plot.title=element_text(hjust=0.5),
axis.title.x=element_text(face="bold",
color="Black",size=12),
axis.title.y=element_text(face="bold",
color="Black",
size=12),
legend.title=element_text(face="bold",
size=10))
Post hoc Tukey
summary(glht(cox_fullmodel4,mcp(Treatment="Tukey")))
Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: Tukey Contrasts
Fit: coxme(formula = Surv(NumDaysAlive, CensorReproduction) ~ Treatment +
(1 | Site/Population) + (1 | Block), data = mydata)
Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
Biomass Removal - Grassland == 0 1.58171 0.22789 6.941 < 0.001 ***
Hemizonia - Grassland == 0 -0.03584 0.24974 -0.144 0.99990
Raking - Grassland == 0 -0.27778 0.30542 -0.909 0.89129
Raking & Clipping - Grassland == 0 0.94460 0.24661 3.830 0.00117 **
Hemizonia - Biomass Removal == 0 -1.61755 0.23484 -6.888 < 0.001 ***
Raking - Biomass Removal == 0 -1.85949 0.28678 -6.484 < 0.001 ***
Raking & Clipping - Biomass Removal == 0 -0.63711 0.21317 -2.989 0.02288 *
Raking - Hemizonia == 0 -0.24193 0.31280 -0.773 0.93699
Raking & Clipping - Hemizonia == 0 0.98044 0.25237 3.885 < 0.001 ***
Raking & Clipping - Raking == 0 1.22237 0.29585 4.132 < 0.001 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)
Plot by 2 treatments (Control and Biomass Removal), and adding in Habitat
cox_sub_mydata<-subset(mydata,Treatment%in%c('Control','Biomass Removal'))
cox_model_graph3<-survfit(Surv(NumDaysAlive,CensorReproduction)~
Treatment,
data=cox_sub_mydata)
autoplot(cox_model_graph3)+
labs(x="\n Survival Time (Days)",
y="Survival Probabilities\n",
title="Survival Times Of \n Roadside and Vegetated Populations\n")+
theme(plot.title=element_text(hjust=0.5),
axis.title.x=element_text(face="bold",
color="Black",
size=12),
axis.title.y=element_text(face="bold",
color="Black",
size=12),
legend.title=element_text(face="bold",
size=10))
cox_model_graph3<-survfit(Surv(NumDaysAlive,CensorReproduction)~
Treatment+
Habitat,
data=cox_sub_mydata)
autoplot(cox_model_graph3)+
labs(x="\n Survival Time (Days)",
y="Survival Probabilities\n",
title="Survival Times Of \n Roadside and Vegetated Populations\n")+
theme(plot.title=element_text(hjust=0.5),
axis.title.x=element_text(face="bold",
color="Black",
size=12),
axis.title.y=element_text(face="bold",
color="Black",
size=12),
legend.title=element_text(face="bold",
size=10))